C C ************************************************************* C * * C * ********************* * C * ********************* * C * ** ** * C * **PROGRAM GEBAT-V ** * C * ** ** * C * ********************* * C * ********************* * C * * C * GENERAL BUNDLE ADJUSTMENT TRIANGULATION * C * * C * * C * BY : SABRY EL-HAKIM..NRC1980 * C * ---------------------------- * C * * C * * C * * C * PROGRAM FEATURES * C * ---------------- * C * -BUNDLE ADJUSTMENT WITH SELF CALIBRATION USING * C * HARMONIC FUNCTION * C * -BLOCK - VARIANT APPROACH * C * -COMBINED ADJUSTMENT WITH SPATIAL DISTANCES * C * AND HEIGHT DIFFERENCES * C * -GROSS-ERROR DETECTION BY DATA SNOOPING * C * -ERROR ELLIPSOIDS ARE COMPUTED * C * * C * LIMITATIONS : ( CAN BE CHANGED BY THE USER ) * C * ------------- * C * -MAX. NO. OF PHOTOS.....................8 * C * -MAX. NO. OF IMAGE PTS................800 * C * -MAX. NO. OF OBJECT PTS...............250 * C * -MAX. NO. OF GEODETIC PTS..............80 * C * -UNLIMITED ; PTS./PHOTO ,CONTROL PTS. * C * & NO. OF DISTANCES&HEIGHTS * C * * C ************************************************************* C C IMPLICIT REAL*8(A-H,O-Z) DIMENSION L(801),LG(801),NF(250),NV(250),LC(250),LQ(250), *ICAT(102),INOT(801) REAL*4 XP(801),YP(801),PX1(8,17),PX2(1,750) *,X1(80,80),X2(80,80),X3(80,80),X4(80,80) REAL*8 XG(801),YG(801),ZG(801),XC(250),YC(250),ZC(250),XQ(250), *YQ(250),ZQ(250),X(136),B(136),P(750) *,UG(240),XO(240),NOG(28920),AG(1,240) REAL*8 AP1(2,17),MP(2,2),MPA(2,17),W(2),UP1(136),XP21(750,17), *NP11(17,136),M(3,3),BP(2,2),AP2(2,3),NP21(750,17),NP22(3,750), *UP2(750),NE(136,17),VEC(750),AL(801),R(1,136),C(9316) 3,FF(10),PF(10) COMMON NPHO(8),IMAX(8),IMIN(8) DATA NOB,NCO,NCH,IL1/250,60,180,21/ DATA IRC,IWP,IL2,IUY,INPT,IUNK,ITUN,IPTO/5,6,23,9316,801,136,136, *750/ DATA IGM,IGN,JGN/80,28920,240/ WRITE(IWP,10) 10 FORMAT(1H1) CALL GEBAT1(L,LG,NF,NV,LC,LQ,XP,YP,XG,YG,ZG,PX2,PX1,XC,YC,ZC,XQ,YQ *,ZQ,X,B,P,AP1,MP,MPA,W,UP1,IWP,IL2,IUY,INPT,IUNK,ITUN,IPTO,NOB,NCO 2,NP11,M,BP,AP2,NP21,NP22,UP2,NE,VEC,ICAT,AL,IRC,C,FF,PF,R,NCH, *XP21,INOT,X1,X2,X3,X4,UG,XO,NOG,AG,IGM,IGN,JGN,IL1) STOP END SUBROUTINE GEBAT1(L,LG,NF,NV,LC,LQ,XP,YP,XG,YG,ZG,PX2,PX1,XC,YC,ZC *,XQ,YQ,ZQ,X,B,P,AP1,MP,MPA,W,UP1,IWP,IL2,IUY,INPT,IUNK,ITUN,IPTO, *NOB,NCO,NP11,M,BP,AP2,NP21,NP22,UP2,NE,VEC,ICAT,AL,IRC,C,FF,PF,R, *NCH,XP21,INOT,X1,X2,X3,X4,UG,XO,NOG,AG,IGM,IGN,JGN,IL1) IMPLICIT REAL*8(A-H,O-Z) INTEGER Q,ICAT(ITUN) DIMENSION L(INPT),LG(INPT),NF(NCO),NV(NCO),LC(NOB),LQ(NOB) *,NPF(8),INOT(INPT) REAL*4 XP(INPT),YP(INPT),PX1(8,17),PX2(1,IPTO) *,X1(IGM,IGM),X2(IGM,IGM),X3(IGM,IGM),X4(IGM,IGM) REAL*8 XG(INPT),YG(INPT),ZG(INPT),XC(NOB),YC(NOB),ZC(NOB), 2XQ(NCH),YQ(NCH),ZQ(NCH),X(ITUN),B(ITUN),P(IPTO),AP1(2,17),W(2), *MP(2,2),MPA(2,17),UP1(ITUN),NP11(17,IUNK),M(3,3),XP21(IPTO,17), 5BP(2,2),AP2(2,03),NP21(IPTO,17),NP22(3,IPTO),UP2(IPTO),NE(ITUN,17) REAL*8 VEC(IPTO),AL(INPT),R(1,ITUN),C(IUY),FF(10),PF(10) *,UG(JGN),XO(JGN),AG(1,JGN),NOG(IGN) COMMON NPHO(8),IMAX(8),IMIN(8) C C ************************************************************* C * * C * PROGRAM LIMITATIONS CAN BE CHANGED BY CHANGING THE * C * PROGRAM DIMENSIONS AND THE FOLLOWING DATA STATEMENT * C * PARAMETERS :- * C * * C * -IUY = (NF*17)**2/2 +(NF*17)/2 * C * WHERE NF : #OF PHOTOS , * C * * C * -IGM = # OF GEODETIC PTS. * C * -IGN =((3*IGM)**2+3*IGM)/2 * C * -JGN = 3*IGM * C * * C * -INPT = # OF IMAGE PTS. (TOTAL) + 1 * C * -IUNK = 17*NF * C * -ITUN = IUNK * C * -NOB = # OF OBJECT PTS. * C * -IPTO = 3*NOB * C * -NCO = # OF CONTROL PTS. * C * -NCH = # OF CHECK PTS. * C * * C * NOTE: THE DIMENSIONS FOR THE ARRAYS NPHO IMAX IMIN * C * IN BOTH MAIN AND A FEW SUBROUTINES SHOULD BE CHANGED * C * ACCORDING TO THE CHANGED NUMBER OF PHOTOS. * C * THE ARRAY NPF SHOULD ALSO BE CHANGED IN DIMENSION C * THE RELATED SUBROUTINES ARE GEBAT1,VNORM1,VNORM2, C * UPEC MORMAL. C * * C * FOR THE ABOVE DIMENSIONS, THE PROGRAM NEEDS 160K +2 FILES* C * (SEE DIMENSION STATEMENTS IN SUBROUTINE GEBATM FOR THE * C * LOCATION OF THE ABOVE VARIABLES AS SUBSCRIPTS..) * C * * C * * C * OTHER DATA STATEMENT PARAMETERS :- * C * -IRC : READ UNIT NUMBER * C * -IWP : WRITE UNIT NUMBER * C * -IL1 : ADDITIONAL ON-LINE FILE NUMBER * C * -IL2 : ADDITIONAL ON-LINE FILE NUMBER * C * * C ************************************************************* C C 1 FORMAT(10A8) 12 FORMAT(I10) 30 FORMAT(5X,I6,19X,2F10.3) 1234 FORMAT(/10X,'PHOTO NO.',I6/9X,17('-')/) 39 FORMAT(I15,2F15.3) 409 FORMAT(3F10.4) 47 FORMAT(I10,3F20.9) 489 FORMAT(2I8,F8.3,2I8,F8.3,2I8,F8.3) 201 FORMAT(6F10.4/4F20.9) 410 FORMAT(4I10,F10.1) 419 FORMAT(3F15.10) 3030 FORMAT(//5X,'TOTAL NO. OF PHOTO PTS. =',I4/5X, * 'TOTAL NO. OF PHOTOS =',I4////5X, *'THE CONTROL POINTS COORDS. :-'/5X,29('-')/) 3103 FORMAT(/5X,'****ERROR****POINT',I6,'HAS NO APPROX. COORDS.****') ICOX=0 NPPH=0 READ(IRC,410)IG,IQ,IGCO,KG,FSCALE READ(IRC,410)NITR,Q,MBW,IVCEE 2345 FORMAT(5X,16('*')/5X,'*',14X,'*'/5X,'*',' READ-IN DATA ','*'/5X, *'*',14X,'*'/5X,16('*')/) WRITE(IWP,2345) IF(IG.EQ.0)GO TO 1991 READ(IRC,2121)(LQ(KK),KK=1,IG) 1991 CONTINUE READ(IRC,1)(PF(I),I=1,10) READ(IRC,1)(FF(I),I=1,10) IB=IG IK=0 J=1 I=1 16 K=1 C READ(IRC,PF)NFIII,XPP,YPP,IRX,IRY IF(IRX.EQ.0)IRX=1 IF(IRY.EQ.0)IRY=1 WRITE(IWP,1234)NFIII NPHO(J)=NFIII NFIII=J 14 READ(IRC,FF)L(I),XP(I),YP(I) IF(L(I).EQ.666666)GO TO 13 IF(L(I).EQ.-66)GO TO 15 WRITE(IWP,39)L(I),XP(I),YP(I) XP(I)=XP(I)-XPP YP(I)=YP(I)-YPP XP(I)=XP(I)/1000.*IRX YP(I)=YP(I)/1000.*IRY LOIO=L(I) IF(IG.EQ.0)GO TO 987 DO 1987 IY=1,IG IF(L(I).EQ.LQ(IY))GO TO 985 1987 CONTINUE 987 L(I)=L(I)+NFIII*1000000 III=0 IF(I.EQ.1)GO TO 190 I2=I-1 DO 913 II=1,I2 LL=(AL(II)+0.01)/1000000000 A1L=LL AR=AL(II)-A1L*1000000000.0+0.01 IR=AR LOIIO=IR/1000 IR=IR-IR/1000*1000 IF(LOIO.EQ.LOIIO)IM=IR IF(LOIO.EQ.LOIIO)III=1 913 CONTINUE 190 IF(III.EQ.0)IB=IB+1 IF(III.EQ.0)IM=IB AL(I)=L(I) AL(I)=AL(I)*1000+IM INOT(I)=IM GO TO 986 985 CONTINUE L(I)=L(I)+NFIII*1000000 III=0 IF(I.EQ.1)GO TO 19 I2=I-1 DO 113 II=1,I2 LL=(AL(II)+0.01)/1000000000 A1L=LL AR=AL(II)-A1L*1000000000.0+0.01 IR=AR LOIIO=IR/1000 IR=IR-IR/1000*1000 IF(LOIO.EQ.LOIIO)IM=IR IF(LOIO.EQ.LOIIO)III=1 113 CONTINUE 19 IF(III.EQ.0)IK=IK+1 IF(III.EQ.0)IM=IK AL(I)=L(I) AL(I)=AL(I)*1000+IM 986 CONTINUE INOT(I)=IM I=I+1 K=K+1 GO TO 14 13 J=J+1 NPF(J-1)=I-1-NPPH NPPH=NPF(J-1)+NPPH GO TO 16 15 NP=I-1 NPF(J)=NP-NPPH C C NP : THE TOTAL NO. OF PHOTO POINTS C NFI=J C C NFI : THE NUMBER OF PHOTOS C WRITE(IWP,3030)NP,NFI DO 1550 I=1,NP MF=L(I)/1000000 L1=L(I)-MF*1000000 DO 1551 J=1,NP L2=L(J)-L(J)/1000000*1000000 IF(I.EQ.J)GO TO 1551 IF(L1.EQ.L2)GO TO 1550 IF(J.EQ.NP)WRITE(IWP,1552)L1,NPHO(MF) IF(J.EQ.NP)ICOX=1 1551 CONTINUE 1550 CONTINUE 1552 FORMAT(/5X,'****ERROR****POINT',I7,' IS SINGLE , REMOVE FROM PHOTO * NO.',I7) DO 6266 J=1,NFI IMIN(J)=10000 6266 IMAX(J)=1 DO 6166 I=1,NP IX=(AL(I)+0.01)/1000000000 AA=IX AR=AL(I)-AA*1000000000.+0.01 ISO=AR KN=ISO-ISO/1000*1000 IF(KN.LT.IMIN(IX))IMIN(IX)=KN IF(KN.GT.IMAX(IX))IMAX(IX)=KN 6166 CONTINUE DO 6366 I=1,NFI IMIN(I)=(IMIN(I)-1)*3+1 6366 IMAX(I)=IMAX(I)*3 READ(IRC,1)(PF(I),I=1,10) I=0 55 I=I+1 IJ=3*I READ(IRC,PF)LC(I),XC(I),YC(I),ZC(I),UP2(IJ-2),UP2(IJ-1),UP2(IJ) IF(LC(I).NE.-66)GO TO 55 IC=I-1 ICR=IC 56 FORMAT(I10) DO 752 I=1,IC 752 LQ(I)=0 DO 762J=1,NP LL=(AL(J)+0.01)/1000000000 A1L=LL AR=AL(J)-A1L*1000000000.0+0.01 JR=AR IR=JR/1000 IN=JR-IR*1000 DO 31 I=1,IC IF(IR.EQ.LC(I))GO TO 32 GO TO 3102 32 XG(J)=XC(I) YG(J)=YC(I) ZG(J)=ZC(I) IJ=3*I PX2(1,3*IN-2)=UP2(IJ-2) PX2(1,3*IN-1)=UP2(IJ-1) PX2(1,3*IN)=UP2(IJ) LQ(I)=LQ(I)+1 GO TO 762 3102 IF(I.EQ.IC)XG(J)=-99999999.9 31 CONTINUE IF(I.EQ.IC)LC(I+1)=IR IF(I.EQ.IC)XC(I+1)=XG(J) IF(I.EQ.IC)IC=IC+1 762 CONTINUE C C READING THE WT. MATRIX OF THE PHOTO POINTS & PX1, PX1U&PX2 C READ(IRC,24)P1,P2,P3 24 FORMAT(3F15.10) PX=P2/(P1*P2-P3**2) PY=P1/(P1*P2-P3**2) PXY=P3/(P1*P2-P3**2) JC=IC NU=17 N1=0 DO 710 II=1,NFI READ(IRC,201)(PX1(II,I),I=12,17),(PX1(II,J),J=1,3),P8 DO 1106 IK=4,11 1106 PX1(II,IK)=100.0 710 CONTINUE READ(IRC,419)PXCO,PYCO,PZCO NUN=NU*NFI NAN=NUN+N1 DO 21 I=1,NAN B(I)=0.0 21 X(I)=0.0 DO 98 I=1,NFI READ(IRC,47)ND,XX,YY,ZZ READ(IRC,489)IOM,JOM,SOM,IPHI,JPHI,SPHI,IKP,JKP,SKP CALL DMSRAD(IOM,JOM,SOM,AOM) CALL DMSRAD(IPHI,JPHI,SPHI,PHI) CALL DMSRAD(IKP,JKP,SKP,AKP) READ(IRC,198)XO2,YO2,F2 NUII=NU*I XO2=XO2/1000. YO2=YO2/1000. F2=F2/1000. X(NUII-16)=XO2 X(NUII-15)=YO2 X(NUII-14)=F2 X(NUII-2)=AOM X(NUII-1)=PHI X(NUII)=AKP X(NUII-5)=XX X(NUII-4)=YY 98 X(NUII-3)=ZZ 198 FORMAT(3F10.3) ITR=0 JCS=JC*3 DO 3698 J=1,JCS 3698 VEC(J)=0.0 READ(IRC,56)ICON READ(IRC,2121)(NF(I),I=1,ICON) 2121 FORMAT(8I10/8I10/8I10/8I10/8I10/8I10/8I10/8I10/8I10/8I10/8I10) DO 6789 I=1,ICON INF=NF(I) DO 982 J=1,IC 982 IF(LC(J).EQ.NF(I))INP=J 6789 WRITE(IWP,39)INF,XC(INP),YC(INP) READ(IRC,56)IVER READ(IRC,2121)(NV(I),I=1,IVER) DO 6790 I=1,IVER INV=NV(I) DO 981 J=1,IC 981 IF(LC(J).EQ.NV(I))INF=J 6790 WRITE(IWP,4747)INV,ZC(INF) 4747 FORMAT(I15,30X,F15.3) WRITE(IWP,782) 782 FORMAT(//5X,'APPROXIMATE OBJECT COORDS.'/5X,26('-')//11X,'PT.', *11X,'XG',13X,'YG',13X,'ZG',6X,'# OF IMAGES'/) 722 FORMAT(I15,3F15.3,I10) DO 742 I=1,ICR WRITE(IWP,722)LC(I),XC(I),YC(I),ZC(I),LQ(I) IF(LQ(I).EQ.0)WRITE(IWP,772)LC(I) IF(LQ(I).EQ.0)ICOX=1 742 CONTINUE 772 FORMAT(/5X,'****ERROR****WRONG PT.(',I6,') IN OBJECT COORDS.LIST') IF(ICOX.EQ.1)STOP I=0 879 I=I+1 READ(IRC,PF)LQ(I),XQ(I),YQ(I),ZQ(I) IF(LQ(I).NE.-66)GO TO 879 ICC=I-1 1531 CONTINUE DO 1122 I=1,ICON DO 980 J=1,NP LL=(AL(J)+0.01)/1000000000 A1L=LL AR=AL(J)-A1L*1000000000.0+0.01 JR=AR IR=JR/1000 980 IF(IR.EQ.NF(I))INF=JR-IR*1000 PX2(1,3*INF-2)=PXCO 1122 PX2(1,3*INF-1)=PYCO DO 3344 I=1,IVER DO 983 J=1,NP LL=(AL(J)+0.01)/1000000000 A1L=LL AR=AL(J)-A1L*1000000000.0+0.01 JR=AR IR=JR/1000 983 IF(IR.EQ.NV(I))INV=JR-IR*1000 3344 PX2(1,3*INV)=PZCO 501 CONTINUE DO 142 J=1,NUN UP1(J)=0.0 DO 142 I=1,NU 142 NP11(I,J)=0.0 9182 CONTINUE ITR=ITR+1 IF(ITR.EQ.Q)GO TO 7614 GO TO 7514 7614 DO 7414 IN=4,11 DO 7414 II=1,NFI 7414 PX1(II,IN)=P8 7514 CONTINUE INB=NP NPS=INB NP=INPT INU=NUN INC=JCS JCS=IPTO NUN=IUNK NAN=ITUN C IF(IG.EQ.0.OR.ITR.GT.NITR)GO TO 228 IF(IGCO.EQ.2)GO TO 328 CALL SNORM(NOG,UG,NAN,NP,XG,YG,ZG,IG,ITR,X1,X2,INB,IGM,L,XO,IGN, *AG,WG,AL,JGN,IRC,IWP) IF(IGCO.EQ.1)GO TO 228 328 CONTINUE CALL DNORM(NOG,UG,NAN,NP,XG,YG,ZG,IG,ITR,X3,X4,INB,IGM,L,XO,IGN, *AG,WG,AL,JGN,IRC,IWP,IGCO) 228 CONTINUE C CALL VNORM1(L,XP,YP,XG,YG,ZG,X,NU,NP,NFI,AP1,MP,W,AL,IWP,IL2, *NP11,INB,MPA,N1,NUN,NAN,UP1,LG,M,BP,PX,PY,PXY,AP2,NP21,JCS,IRC, *NP22,UP2,MBW,B,VEC,ITR,INC,INU,PX1,F2,NITR,C,IUY,CEGX,CEGY,FSCALE) IF(ITR.LE.NITR)GO TO 1956 IF(IVCEE.EQ.0)RETURN CALL VNORM2(L,XP,YP,XG,YG,ZG,X,NU,NP,NFI,AP1,MP,W,AL,IWP,IL2, *NP11,INB,MPA,N1,NUN,NAN,UP1,LG,M,BP,PX,PY,PXY,AP2,NP21,JCS,IRC, *NP22,UP2,MBW,B,VEC,ITR,INC,INU,PX1,F2,NITR,C,IUY,CEGX,CEGY,ICON, *IVER) RETURN 1956 CONTINUE NP=INB CALL NORM2(NP22,NOG,PX2,JCS,JC,IG,P,M,NAN,3,INC,INU,IGN,IWP,IL1,R) 7723 CONTINUE CALL UPEC(UP1,INU,NP21,UP2,NP22,P,JC,NAN,JCS,VEC,3,IL2,INC,NOG, *UG,IG,IGN,JGN) 7724 CONTINUE ICR=0 DO 7714 I=1,NFI DO 7714 J=1,NU ICR=ICR+1 7714 UP1(ICR)=UP1(ICR)-B(ICR)*PX1(I,J) CALL MORMAL(JC,JCS,NUN,NAN,NP22,NP11,PX1,NE,VEC,NP21,INC,INU, *IL2,C,IUY,IWP,ICAT,XP21,ITR,NITR,NPF,INOT,INPT,IG,NOG,IGN,NU) NAM=NAN JCS=INC NUN=INU NAN=NUN+N1 CALL QHLSK(NAM,C,R,INDEF,ITR,N1,INU,Q,IUY,IWP,ICAT) IF(ITR.LT.NITR)GO TO 4260 CALL UTINV(C,NAM,INU,P,IUY) REWINDIL1 DO 10 I=1,INU DO 11 K=I,INU P(K)=0.0 DO 11 J=K,INU IOJ=((J-1)**2+J-1)/2+I KOJ=((J-1)**2+J-1)/2+K 11 P(K)=P(K)+C(IOJ)*C(KOJ) WRITE(IL1)(P(II),II=I,INU) 10 CONTINUE ENDFILEIL1 REWINDIL1 DO 90 I=1,INU READ(IL1)(P(II),II=I,INU) DO 90 J=I,INU IOJ=((J-1)**2+J-1)/2+I 90 C(IOJ)=P(J) GO TO 501 4260 CALL QHLSKS(NAM,C,UP1,B,P,ITR,N1,INU,Q,IUY,ICAT) DO 427 I=1,NAN 427 X(I)=X(I)+B(I) 22 CONTINUE DO 320 I=1,NFI M1=I*17-5 M2=M1-1+2 M3=M2-1+2 M4=M3+1 M5=M4+1 M6=M5+1 WRITE(IWP,330) DEG1=X(M4)*18.D1/3.141592653589793D0 DEG2=X(M5)*18.D1/3.141592653589793D0 DEG3=X(M6)*18.D1/3.141592653589793D0 DEG4=B(M4)*18.D1/3.141592653589793D0 DEG5=B(M5)*18.D1/3.141592653589793D0 DEG6=B(M6)*18.D1/3.141592653589793D0 WRITE(IWP,340)K,X(M1),X(M2),X(M3),DEG1,DEG2,DEG3 WRITE(IWP,341)B(M1),B(M2),B(M3),DEG4,DEG5,DEG6 WRITE(IWP,342) K=I*NU-17 345 WRITE(IWP,350)(X(K+J),B(K+J),J=1,11) 320 CONTINUE 350 FORMAT(32X,'CORR.'/5X,'XO =',2F15.6/5X,'YO =',2F15.6/5X,'F =', *2F15.6/5X,'A00=',2F15.6/5X,'A11=',2F15.6/5X,'B11=',2F15.6/5X,'A20= *',2F15.6/5X,'A22=',2F15.6/5X,'B22=',2F15.6/5X,'A31=',2F15.6/5X, *'B31=',2F15.6) 340 FORMAT(I8,'-ADJ.PAR.=',3F12.3,3F12.4) 341 FORMAT(9X,'CORRECN.=',3F12.3,3F12.4//) 330 FORMAT(25X,'XC',10X,'YC',10X,'ZC',9X,'OMEGA',9X,'PHI',7X,'KAPPA'/) 342 FORMAT(5X,'THE BASIC PAR. + HARMONIC PAR.'/5X,32('-')) REWINDIL2 DO 40 I=1,JCS 40 P(I)=0.0 IK=0 DO 402 J=1,NFI I1=IMIN(J) I2=IMAX(J) DO 402 I=1,NU IK=IK+1 READ(IL2)(NP21(KK,I),KK=I1,I2) DO 402 K=I1,I2 402 P(K)=P(K)+NP21(K,I)*B(IK) DO 709 I=1,JCS 709 UP2(I)=UP2(I)+P(I) JCS=IPTO NUN=IUNK NAN=ITUN CALL UPEC(UP1,INU,NP21,UP2,NP22,P,JC,NAN,JCS,VEC,4,IL2,INC,NOG, *UG,IG,IGN,JGN) 1072 CONTINUE JCS=INC NUN=INU NAN=NUN+N1 WRITE(IWP,20) 20 FORMAT(///5X,'ADJUSTED COORDS. & CORRECTIONS'/5X,30('*')) 3330 FORMAT(/I10,'-',3F12.3/11X,3F12.3) DO 780 I=1,NP LL=(AL(I)+0.01)/1000000000 A1L=LL AR=AL(I)-A1L*1000000000.0+0.01 KN=AR KN=KN-KN/1000*1000 LGOIO=AR/1000 XG(I)=XG(I)-VEC(3*KN-2) YG(I)=YG(I)-VEC(3*KN-1) ZG(I)=ZG(I)-VEC(3*KN) IF(I.EQ.1)GO TO 1003 I2=I-1 DO 1004 JJ=1,I2 LL=(AL(I)+0.01)/1000000000 A1L=LL AR=AL(I)-A1L*1000000000.0+0.01 LGOIO=AR LGOIO=LGOIO/1000 LL=(AL(JJ)+0.01)/1000000000 A1L=LL AR=AL(JJ)-A1L*1000000000.0+0.01 LGOJJO=AR LGOJJO=LGOJJO/1000 IF(LGOIO.EQ.LGOJJO)GO TO 780 1004 CONTINUE 1003 WRITE(IWP,3330)LGOIO,XG(I),YG(I),ZG(I),VEC(3*KN-2), 1VEC(3*KN-1),VEC(3*KN) 780 CONTINUE IF(IG.EQ.0)GO TO 881 DO 882 I=1,IG I3=I*3 XO(I3)=XO(I3)-VEC(I3) XO(I3-1)=XO(I3-1)-VEC(I3-1) 882 XO(I3-2)=XO(I3-2)-VEC(I3-2) 881 CONTINUE NP=INPT CALL CHECK(LG,XG,YG,ZG,LQ,XQ,YQ,ZQ,NQ,NP,INB,ICC,AL,NCH,IWP) NP=INB 5011 IF(ITR.LT.NITR)GO TO 501 IF(IG.EQ.0)STOP IF(IGCO.EQ.2)GO TO 428 CALL SNORM(NOG,UG,NAN,NP,XG,YG,ZG,IG,0,X1,X2,INB,IGM,L,XO,IGN,AG, *WG,AL,JGN,IRC,IWP) IF(IGCO.EQ.1)RETURN 428 CONTINUE CALL DNORM(NOG,UG,NAN,NP,XG,YG,ZG,IG,0,X3,X4,INB,IGM,L,XO,IGN,AG, *WG,AL,JGN,IRC,IWP,IGCO) RETURN END SUBROUTINE VNORM1(L,XP,YP,XG,YG,ZG,X,NU,NQ,NFI,AP1,MP,W,AL,IWP,IL2 *,NP11,INB,MPA,N1,NUM,NAM,UP1,LG,M,BP,PX,PY,PXY,AP2,NP21,JJS,IRC, *NP22,UP2,MBW,B,VEC,ITR,INC,INU,PX1,F2,NITR,C,MOM,CEGX,CEGY,FSCALE) C C IMPLICIT REAL*8(A-H,O-Z) DIMENSION L(NQ),XG(NQ),YG(NQ),ZG(NQ),LG(NQ) REAL*4 XP(NQ),YP(NQ),PX1(8,17) REAL*8 AP1(2,17),MP(2,2),MPA(2,17),W(2),UP1(NAM),NP22(3,JJS), *NP11(17,NUM),M(3,3),BP(2,2),NP21(JJS,17),VEC(JJS),AP2(2,3),AL(NQ), *B(NAM),X(NAM),UP2(JJS),C(MOM) DATA PASS,FAIL,UNDIF/'O.K.','REJECT','U.D.??'/ COMMON NPHO(8),IMAX(8),IMIN(8) IZO=0 WRITE(IWP,51)ITR 50 FORMAT(//28X,'RESIDUALS',8X,'REDUNDANCY NOS.', *5X,'DATA SNOOPING TEST(C=',F4.2,')'/ *28X,9('-'),8X,15('-'),5X,18('-')/11X,'PHOTO-PT.',7X,'VX',8X,'VY', *7X,'RX',8X,'RY',8X,'WX',8X,'WY'/) 51 FORMAT(1H1//5X,'*RESULTS OF ITERATION #',I2/5X,25('*')/) 2458 FORMAT(/) 2973 FORMAT(7X,I6,'-',I6,2F10.1,4F10.3,' <---',A8,F10.3,F8.3,'*',A8) 2468 FORMAT(F10.3) IF(ITR.EQ.1)READ(IRC,2468)TOL IF(MBW.GE.0.AND.ITR.GT.NITR)WRITE(IWP,50)TOL QVX=0.0 QVY=0.0 CHX=0.0 CHY=0.0 JCS=INC NUN=INU NAN=NUN NP=INB REWINDIL2 INDX=1 SEGX=0.0 SEGY=0.0 DO 365 I=1,JCS DO 365 J=1,NU 365 NP21(I,J)=0.0 IF(ITR.GT.NITR)GO TO 200 DO1262 J=1,JCS 1262 UP2(J)=0.0 DO 262 J=1,JCS DO 262 I=1,3 262 NP22(I,J)=0.0 200 J=1 SO=DSIN(X(15)) CO=DCOS(X(15)) SF=DSIN(X(16)) CF=DCOS(X(16)) SK=DSIN(X(17)) CK=DCOS(X(17)) M(1,1)=CF*CK M(1,2)=CO*SK+SO*SF*CK M(1,3)=SO*SK-SF*CK*CO M(2,1)=-CF*SK M(2,2)=CO*CK-SO*SF*SK M(2,3)=SO*CK+SF*SK*CO M(3,1)=SF M(3,2)=-SO*CF M(3,3)=CO*CF 100 IX=(AL(J)+0.01)/1000000000 TEST=PASS NUONFI=NUN NUOIX=NU*IX IIO=0 IF(IX.GT.INDX)IIO=1 IF(IX.GT.INDX)INDX=IX XR=XP(J)-X(NUOIX-16) YR=YP(J)-X(NUOIX-15) R=DSQRT(XR**2+YR**2) KO=IX*NU IF(IIO.EQ.0)GO TO 101 SO=DSIN(X(NUOIX-2)) CO=DCOS(X(NUOIX-2)) CF=DCOS(X(NUOIX-1)) SF=DSIN(X(NUOIX-1)) CK=DCOS(X(NUOIX)) SK=DSIN(X(NUOIX)) M(1,1)=CF*CK M(1,2)=CO*SK+SO*SF*CK M(1,3)=SO*SK-SF*CK*CO M(2,1)=-CF*SK M(2,2)=CO*CK-SO*SF*SK M(2,3)=SO*CK+SF*SK*CO M(3,1)=SF M(3,2)=-SO*CF M(3,3)=CO*CF 101 CONTINUE IF(XG(J).LT.-99999998.)GO TO 401 GO TO 402 401 IF(IX.EQ.1)GO TO 405 IJ=J-1 AF=AL(J)-IX*1000000000.+.01 IA=AF DO 404 ICX=1,IJ IIX=(AL(ICX)+0.01)/1000000000 AB=AL(ICX)-IIX*1000000000.0+0.01 IB=AB IF(IA.EQ.IB)GO TO 406 404 CONTINUE 405 CONTINUE NUNX3=NUOIX-14 XG(J)=X(NUOIX-5)+FSCALE*(M(1,1)*XR+M(1,2)*YR-X(NUNX3)*M(1,3)) YG(J)=X(NUOIX-4)+FSCALE*(M(2,1)*XR+M(2,2)*YR-X(NUNX3)*M(2,3)) ZG(J)=X(NUOIX-3)+FSCALE*(M(3,1)*XR+M(3,2)*YR-X(NUNX3)*M(3,3)) GO TO 402 406 XG(J)=XG(ICX) YG(J)=YG(ICX) ZG(J)=ZG(ICX) 402 CONTINUE T1=M(1,1)*(XG(J)-X(NUOIX-5))+M(1,2)*(YG(J)-X(NUOIX-4))+M(1,3)* *(ZG(J)-X(NUOIX-3)) T2=M(2,1)*(XG(J)-X(NUOIX-5))+M(2,2)*(YG(J)-X(NUOIX-4))+M(2,3)* *(ZG(J)-X(NUOIX-3)) T3=M(3,1)*(XG(J)-X(NUOIX-5))+M(3,2)*(YG(J)-X(NUOIX-4))+M(3,3)* *(ZG(J)-X(NUOIX-3)) 1110 SL=YR/R CL=XR/R S2L=2.*SL*CL C2L=2.*CL**2-1.0 NU=NUOIX-17 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* *S2L+X(NU+10)*R**2*CL+X(NU+11)*R**2*SL NK=-3 NU=17 DVX=T*XR DVY=T*YR KO=NU NUONFI=NUOIX-17 AP1(1,KO-5)=-(XR+DVX)*M(3,1)-X(NUONFI+3)*M(1,1) AP1(2,KO-5)=-(YR+DVY)*M(3,1)-X(NUONFI+3)*M(2,1) AP1(1,KO-4)=-(XR+DVX)*M(3,2)-X(NUONFI+3)*M(1,2) AP1(2,KO-4)=-(YR+DVY)*M(3,2)-X(NUONFI+3)*M(2,2) AP1(1,KO-3)=-(XR+DVX)*M(3,3)-X(NUONFI+3)*M(1,3) AP1(2,KO-3)=-(YR+DVY)*M(3,3)-X(NUONFI+3)*M(2,3) AP1(1,KO-2)=(XR+DVX)*(-M(3,3)*(YG(J)-X(NUOIX-4))+M(3,2)*(ZG(J)- *X(NUOIX-3)))+X(NUONFI+3)*(-M(1,3)*(YG(J)-X(NUOIX-4))+M(1,2)*(ZG(J) *-X(NUOIX-3))) AP1(2,KO-2)=(YR+DVY)*(-M(3,3)*(YG(J)-X(NUOIX-4))+M(3,2)*(ZG(J)- *X(NUOIX-3)))+X(NUONFI+3)*(-M(2,3)*(YG(J)-X(NUOIX-4))+M(2,2)*(ZG(J) *-X(NUOIX-3))) AP1(1,KO-1)=(XR+DVX)*(CF*(XG(J)-X(NUOIX-5))+SO*SF*(YG(J)-X(NUOIX-4 *))-CO*SF*(ZG(J)-X(NUOIX-3)))+X(NUONFI+3)*(-SF*CK*(XG(J)-X(NUOIX-5) *)+SO*CF*(YG(J)-X(NUOIX-4))*CK-CF*CK*CO*(ZG(J)-X(NUOIX-3))) AP1(2,KO-1)=(YR+DVY)*(CF*(XG(J)-X(NUOIX-5))+SO*SF*(YG(J)-X(NUOIX-4 *))-CO*SF*(ZG(J)-X(NUOIX-3)))+X(NUONFI+3)*(SF*SK*(XG(J)-X(NUOIX-5)) *-SO*CF*SK*(YG(J)-X(NUOIX-4))+CF*SK*CO*(ZG(J)-X(NUOIX-3))) AP1(1,KO)=X(NUONFI+3)*(M(2,1)*(XG(J)-X(NUOIX-5))+M(2,2)*(YG(J)- *X(NUOIX-4))+M(2,3)*(ZG(J)-X(NUOIX-3))) AP1(2,KO)=X(NUONFI+3)*(-M(1,1)*(XG(J)-X(NUOIX-5))-M(1,2)*(YG(J)- *X(NUOIX-4))-M(1,3)*(ZG(J)-X(NUOIX-3))) AP1(1,NK+4)=-T3-T*T3 AP1(2,NK+4)=0.0 AP1(1,NK+5)=0.0 AP1(2,NK+5)=AP1(1,NK+4) AP1(1,NK+6)=T1 AP1(2,NK+6)=T2 AP1(1,NK+7)=XR*T3 AP1(2,NK+7)=YR*T3 AP1(1,NK+8)=XR*T3*CL AP1(2,NK+8)=YR*T3*CL AP1(1,NK+9)=XR*T3*SL AP1(2,NK+9)=YR*T3*SL AP1(1,NK+10)=XR*T3*R AP1(2,NK+10)=YR*T3*R AP1(1,NK+11)=XR*T3*R*C2L AP1(2,NK+11)=YR*T3*R*C2L AP1(1,NK+12)=XR*T3*R*S2L AP1(2,NK+12)=YR*T3*R*S2L AP1(1,NK+13)=XR*T3*R**2*CL AP1(2,NK+13)=YR*T3*R**2*CL AP1(1,NK+14)=XR*T3*R**2*SL AP1(2,NK+14)=YR*T3*R**2*SL NUONFI=NUN 300 LL=(AL(J)+0.01)/1000000000 A1L=LL AR=AL(J)-A1L*1000000000.0+0.01 ISO=AR KN=ISO ISO=ISO/1000 KN=KN-KN/1000*1000 IIKN=3 NUNO3=NUOIX-14 AP2(1,IIKN-2)=(XR+DVX)*M(3,1)+X(NUNO3)*M(1,1) AP2(2,IIKN-2)=(YR+DVY)*M(3,1)+X(NUNO3)*M(2,1) AP2(1,IIKN-1)=(XR+DVX)*M(3,2)+X(NUNO3)*M(1,2) AP2(2,IIKN-1)=(YR+DVY)*M(3,2)+X(NUNO3)*M(2,2) AP2(1,IIKN)=(XR+DVX)*M(3,3)+X(NUNO3)*M(1,3) AP2(2,IIKN)=(YR+DVY)*M(3,3)+X(NUNO3)*M(2,3) W1=(XR+DVX)*T3+X(NUNO3)*T1 W2=(YR+DVY)*T3+X(NUNO3)*T2 CX=(W1/T3)*1000000. CY=(W2/T3)*1000000. NON=NUN NUN=NUOIX-17 T5=X(NUN+5)*SL**2/R-X(NUN+6)*CL*SL/R+X(NUN+7)*CL+X(NUN+8)*(CL*C2L+ *2.*S2L*SL)+X(NUN+9)*(CL*S2L-2.*SL*C2L)+X(NUN+10)*(XR*CL+ *R)+X(NUN+11)*XR*SL T6=-X(NUN+5)*SL*CL/R+X(NUN+6)*CL**2/R+X(NUN+7)*SL+X(NUN+8)*(SL*C2L *-2.*S2L*CL)+X(NUN+9)*(SL*S2L+2.*C2L*CL)+X(NUN+10)*XR*SL+ *X(NUN+11)*(R+YR*SL) NUN=NON BP(1,1)=T3+T3*(T+XR*T5) BP(1,2)=T3*(XR*T6) BP(2,1)=T3*YR*T5 BP(2,2)=T3+T3*(T+YR*T6) AMP11=BP(1,1)*(PX*BP(1,1)+PXY*BP(1,2))+BP(1,2)*(PXY*BP(1,1)+PY* *BP(1,2)) AMP12=BP(1,1)*(PX*BP(2,1)+PXY*BP(2,2))+BP(1,2)*(PXY*BP(2,1)+PY* *BP(2,2)) AMP21=BP(2,1)*(PX*BP(1,1)+PXY*BP(1,2))+BP(2,2)*(PXY*BP(1,1)+PY* *BP(1,2)) AMP22=BP(2,1)*(PX*BP(2,1)+PXY*BP(2,2))+BP(2,2)*(PXY*BP(2,1)+PY*BP *(2,2)) DET=AMP11*AMP22-AMP12*AMP21 MP(1,1)=AMP22/DET MP(1,2)=-AMP12/DET MP(2,1)=-AMP21/DET MP(2,2)=AMP11/DET 600 W(1)=W1*MP(1,1)+W2*MP(1,2) W(2)=W1*MP(2,1)+W2*MP(2,2) DO 110 II=1,2 DO 120 K=1,NU MPA(II,K)=0.0 DO 120 JJ=1,2 120 MPA(II,K)=MPA(II,K)+MP(II,JJ)*AP1(JJ,K) 110 CONTINUE IF(IIO.EQ.0)GO TO 5150 I1=IMIN(IX-1) I2=IMAX(IX-1) 1005 DO 6150 KL=1,NU WRITE(IL2)(NP21(JJ,KL),JJ=I1,I2) 6150 CONTINUE IF(IZO.EQ.1)GO TO 1000 DO 363 I=1,JCS DO 363 K=1,NU 363 NP21(I,K)=0.0 5150 CONTINUE DO 150 II=1,3 I2=KN*3-3+II DO 160 K=1,NU DO 161 JJ=1,2 161 NP21(I2,K)=NP21(I2,K)+AP2(JJ,II)*MPA(JJ,K) 160 CONTINUE 150 CONTINUE IF(ITR.LE.NITR)GO TO 6006 CALL COEFS(C,AP1,MP,QVX,QVY,MOM,IX,VEC,JJS,NP22,AP2,KN,NP21) XVX=DABS(QVX) XVY=DABS(QVY) EVX=DSQRT(XVX) EVY=DSQRT(XVY) CHX=CX/(CEGX*EVX) CHY=CY/(CEGY*EVY) OCX=DABS(CHX) OCY=DABS(CHY) IF(OCX.GT.TOL.OR.OCY.GT.TOL)TEST=FAIL IF(QVX.LT.0.0.OR.QVY.LT.0.0)TEST=UNDIF IF(IIO.EQ.1)WRITE(IWP,2458) IF(MBW.NE.0)WRITE(IWP,2973)NPHO(IX),ISO,CX,CY,XVX,XVY,CHX,CHY,TEST GO TO 501 6006 CONTINUE IF(MBW.NE.0)WRITE(IWP,2973)NPHO(IX),ISO,CX,CY DO 130 II=1,NU DO 140 K=1,NU KK=IX*NU-NU+K DO 140 JJ=1,2 140 NP11(II,KK)=NP11(II,KK)+AP1(JJ,II)*MPA(JJ,K) 130 CONTINUE DO 210 II=1,NU IK=IX*NU-NU+II DO 210 KK=1,2 210 UP1(IK)=UP1(IK)+AP1(KK,II)*W(KK) DO 111 II=1,2 DO 121 K=1,3 MPA(II,K)=0.0 DO 121 JJ=1,2 121 MPA(II,K)=MPA(II,K)+MP(II,JJ)*AP2(JJ,K) 111 CONTINUE DO 131 II=1,3 DO 141 K=1,3 KK=KN*3-3+K DO 141 JJ=1,2 141 NP22(II,KK)=NP22(II,KK)+AP2(JJ,II)*MPA(JJ,K) 131 CONTINUE 500 DO 211 II=1,3 IK=KN*3-3+II DO 211 KK=1,2 211 UP2(IK)=UP2(IK)+AP2(KK,II)*W(KK) 501 CONTINUE SEGX=SEGX+CX**2 SEGY=SEGY+CY**2 IF(J.EQ.NP)IZO=1 IF(J.EQ.NP)GO TO 2005 J=J+1 GO TO 100 2005 I1=IMIN(IX) I2=IMAX(IX) GO TO 1005 1000 DF=NP CEGX=DSQRT(SEGX/DF) CEGY=DSQRT(SEGY/DF) WRITE(IWP,1001)CEGX,CEGY 1001 FORMAT(//5X,'ST. ERRORS'/5X,9('-')/5X,'SEGX=',F8.1/5X,'SEGY=', *F8.1//) ENDFILEIL2 IF(ITR.LE.NITR)RETURN RETURN END SUBROUTINE VNORM2(L,XP,YP,XG,YG,ZG,X,NU,NQ,NFI,AP1,MP,W,AL,IWP,IL2 *,NP11,INB,MPA,N1,NUM,NAM,UP1,LG,M,BP,PX,PY,PXY,AP2,NP21,JJS,IRC, *NP22,UP2,MBW,B,VEC,ITR,INC,INU,PX1,F2,NITR,C,MOM,CEGX,CEGY,ICO,IV) C C IMPLICIT REAL*8(A-H,O-Z) DIMENSION L(NQ),XG(NQ),YG(NQ),ZG(NQ),LG(NQ) REAL*4 XP(NQ),YP(NQ),PX1(8,17) REAL*8 AP1(2,17),MP(2,2),MPA(2,17),W(2),UP1(NAM),NP22(3,JJS), *NP11(17,NUM),M(3,3),BP(2,2),NP21(JJS,17),VEC(JJS),AP2(2,3),AL(NQ), *B(NAM),X(NAM),UP2(JJS),C(MOM),AXX(17),V(3,3),ST(3) DATA PASS,FAIL,UNDIF/'O.K.','REJECT','U.D.??'/ DATA WD1,WD2/'STANDARD','NON-STRD'/ COMMON NPHO(8),IMAX(8),IMIN(8) CK=1.0 READ(IRC,1955)ICL IF(ICL.EQ.199)CK=1.0 IF(ICL.EQ.500)CK=1.538 IF(ICL.EQ.900)CK=2.5 IF(ICL.EQ.950)CK=2.7 IF(ICL.EQ.990)CK=3.368 IF(CK.LE.1.0)ICL=199 QK=CK ST1=0.0 ST2=0.0 ST3=0.0 CCL=ICL CCL=CCL/10.0 WD=WD1 IF(CCL.GT.20.0)WD=WD2 SEGX=0.0 SEGY=0.0 SEGZ=0.0 JCS=INC NUN=INU NAN=NUN NP=INB INDX=1 JC=JCS/3 SEUW=(CEGX+CEGY)/2000000 WRITE(IWP,51) WRITE(IWP,1973) II=0 DO 1972 I=1,NFI DO 1974 J=1,NU II=II+1 IOI=((II-1)**2+II-1)/2+II 1974 AXX(J)=DSQRT(C(IOI))*SEUW WRITE(IWP,1975)I,(AXX(JHG),JHG=12,17) 1972 WRITE(IWP,1876)(AXX(KUY),KUY=1,11) 1876 FORMAT(/11X,'XO - YO - F :',3D13.3//10X,'ADDED PARAMETERS :', *D9.3/28X,D9.3/28X,D9.3/28X,D9.3/28X,D9.3/28X,D9.3/28X,D9.3/28X, *D9.3///) 1973 FORMAT(5X,'ST. DIV. OF OR. PARAMETERS'/5X,26('_')/20X,'XC',11X, *'YC',11X,'ZC',11X,'OM',11X,'FI',11X,'KP'/) 1975 FORMAT(I10,'-',6D13.3) II=0 DO 1948 III=1,JC DO 1948 I=1,NP LL=(AL(I)+0.01)/1000000000 A1L=LL AR=AL(I)-A1L*1000000000.0+0.01 ISO=AR KN=ISO-ISO/1000*1000 IF(KN.NE.III)GO TO 1948 II=II+1 LS=L(II) XGS=XG(II) YGS=YG(II) ZGS=ZG(II) ALS=AL(II) XPS=XP(II) YPS=YP(II) L(II)=L(I) XG(II)=XG(I) YG(II)=YG(I) ZG(II)=ZG(I) AL(II)=AL(I) XP(II)=XP(I) YP(II)=YP(I) L(I)=LS XG(I)=XGS YG(I)=YGS ZG(I)=ZGS AL(I)=ALS XP(I)=XPS YP(I)=YPS 1948 CONTINUE IZO=0 2458 FORMAT(/) 2468 FORMAT(F10.3) DO 365 I=1,JCS DO 365 J=1,NU 365 NP21(I,J)=0.0 200 J=1 100 IX=(AL(J)+0.01)/1000000000 NUONFI=NUN NUOIX=NU*IX XR=XP(J)-X(NUOIX-16) YR=YP(J)-X(NUOIX-15) R=DSQRT(XR**2+YR**2) KO=IX*NU SO=DSIN(X(NUOIX-2)) CO=DCOS(X(NUOIX-2)) CF=DCOS(X(NUOIX-1)) SF=DSIN(X(NUOIX-1)) CK=DCOS(X(NUOIX)) SK=DSIN(X(NUOIX)) M(1,1)=CF*CK M(1,2)=CO*SK+SO*SF*CK M(1,3)=SO*SK-SF*CK*CO M(2,1)=-CF*SK M(2,2)=CO*CK-SO*SF*SK M(2,3)=SO*CK+SF*SK*CO M(3,1)=SF M(3,2)=-SO*CF M(3,3)=CO*CF 101 CONTINUE T1=M(1,1)*(XG(J)-X(NUOIX-5))+M(1,2)*(YG(J)-X(NUOIX-4))+M(1,3)* *(ZG(J)-X(NUOIX-3)) T2=M(2,1)*(XG(J)-X(NUOIX-5))+M(2,2)*(YG(J)-X(NUOIX-4))+M(2,3)* *(ZG(J)-X(NUOIX-3)) T3=M(3,1)*(XG(J)-X(NUOIX-5))+M(3,2)*(YG(J)-X(NUOIX-4))+M(3,3)* *(ZG(J)-X(NUOIX-3)) 1110 SL=YR/R CL=XR/R S2L=2.*SL*CL C2L=2.*CL**2-1.0 NU=NUOIX-17 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* *S2L+X(NU+10)*R**2*CL+X(NU+11)*R**2*SL NK=-3 NU=17 DVX=T*XR DVY=T*YR KO=NU NUONFI=NUOIX-17 AP1(1,KO-5)=-(XR+DVX)*M(3,1)-X(NUONFI+3)*M(1,1) AP1(2,KO-5)=-(YR+DVY)*M(3,1)-X(NUONFI+3)*M(2,1) AP1(1,KO-4)=-(XR+DVX)*M(3,2)-X(NUONFI+3)*M(1,2) AP1(2,KO-4)=-(YR+DVY)*M(3,2)-X(NUONFI+3)*M(2,2) AP1(1,KO-3)=-(XR+DVX)*M(3,3)-X(NUONFI+3)*M(1,3) AP1(2,KO-3)=-(YR+DVY)*M(3,3)-X(NUONFI+3)*M(2,3) AP1(1,KO-2)=(XR+DVX)*(-M(3,3)*(YG(J)-X(NUOIX-4))+M(3,2)*(ZG(J)- *X(NUOIX-3)))+X(NUONFI+3)*(-M(1,3)*(YG(J)-X(NUOIX-4))+M(1,2)*(ZG(J) *-X(NUOIX-3))) AP1(2,KO-2)=(YR+DVY)*(-M(3,3)*(YG(J)-X(NUOIX-4))+M(3,2)*(ZG(J)- *X(NUOIX-3)))+X(NUONFI+3)*(-M(2,3)*(YG(J)-X(NUOIX-4))+M(2,2)*(ZG(J) *-X(NUOIX-3))) AP1(1,KO-1)=(XR+DVX)*(CF*(XG(J)-X(NUOIX-5))+SO*SF*(YG(J)-X(NUOIX-4 *))-CO*SF*(ZG(J)-X(NUOIX-3)))+X(NUONFI+3)*(-SF*CK*(XG(J)-X(NUOIX-5) *)+SO*CF*(YG(J)-X(NUOIX-4))*CK-CF*CK*CO*(ZG(J)-X(NUOIX-3))) AP1(2,KO-1)=(YR+DVY)*(CF*(XG(J)-X(NUOIX-5))+SO*SF*(YG(J)-X(NUOIX-4 *))-CO*SF*(ZG(J)-X(NUOIX-3)))+X(NUONFI+3)*(SF*SK*(XG(J)-X(NUOIX-5)) *-SO*CF*SK*(YG(J)-X(NUOIX-4))+CF*SK*CO*(ZG(J)-X(NUOIX-3))) AP1(1,KO)=X(NUONFI+3)*(M(2,1)*(XG(J)-X(NUOIX-5))+M(2,2)*(YG(J)- *X(NUOIX-4))+M(2,3)*(ZG(J)-X(NUOIX-3))) AP1(2,KO)=X(NUONFI+3)*(-M(1,1)*(XG(J)-X(NUOIX-5))-M(1,2)*(YG(J)- *X(NUOIX-4))-M(1,3)*(ZG(J)-X(NUOIX-3))) AP1(1,NK+4)=-T3-T*T3 AP1(2,NK+4)=0.0 AP1(1,NK+5)=0.0 AP1(2,NK+5)=AP1(1,NK+4) AP1(1,NK+6)=T1 AP1(2,NK+6)=T2 AP1(1,NK+7)=XR*T3 AP1(2,NK+7)=YR*T3 AP1(1,NK+8)=XR*T3*CL AP1(2,NK+8)=YR*T3*CL AP1(1,NK+9)=XR*T3*SL AP1(2,NK+9)=YR*T3*SL AP1(1,NK+10)=XR*T3*R AP1(2,NK+10)=YR*T3*R AP1(1,NK+11)=XR*T3*R*C2L AP1(2,NK+11)=YR*T3*R*C2L AP1(1,NK+12)=XR*T3*R*S2L AP1(2,NK+12)=YR*T3*R*S2L AP1(1,NK+13)=XR*T3*R**2*CL AP1(2,NK+13)=YR*T3*R**2*CL AP1(1,NK+14)=XR*T3*R**2*SL AP1(2,NK+14)=YR*T3*R**2*SL NUONFI=NUN 300 LL=(AL(J)+0.01)/1000000000 A1L=LL AR=AL(J)-A1L*1000000000.0+0.01 ISO=AR KN=ISO ISO=ISO/1000 KN=KN-KN/1000*1000 IIO=0 IF(KN.GT.INDX)IIO=1 IF(KN.GT.INDX)INDX=KN IIKN=3 NUNO3=NUOIX-14 AP2(1,IIKN-2)=(XR+DVX)*M(3,1)+X(NUNO3)*M(1,1) AP2(2,IIKN-2)=(YR+DVY)*M(3,1)+X(NUNO3)*M(2,1) AP2(1,IIKN-1)=(XR+DVX)*M(3,2)+X(NUNO3)*M(1,2) AP2(2,IIKN-1)=(YR+DVY)*M(3,2)+X(NUNO3)*M(2,2) AP2(1,IIKN)=(XR+DVX)*M(3,3)+X(NUNO3)*M(1,3) AP2(2,IIKN)=(YR+DVY)*M(3,3)+X(NUNO3)*M(2,3) W1=(XR+DVX)*T3+X(NUNO3)*T1 W2=(YR+DVY)*T3+X(NUNO3)*T2 NON=NUN NUN=NUOIX-17 T5=X(NUN+5)*SL**2/R-X(NUN+6)*CL*SL/R+X(NUN+7)*CL+X(NUN+8)*(CL*C2L+ *2.*S2L*SL)+X(NUN+9)*(CL*S2L-2.*SL*C2L)+X(NUN+10)*(XR*CL+ *R)+X(NUN+11)*XR*SL T6=-X(NUN+5)*SL*CL/R+X(NUN+6)*CL**2/R+X(NUN+7)*SL+X(NUN+8)*(SL*C2L *-2.*S2L*CL)+X(NUN+9)*(SL*S2L+2.*C2L*CL)+X(NUN+10)*XR*SL+ *X(NUN+11)*(R+YR*SL) NUN=NON BP(1,1)=T3+T3*(T+XR*T5) BP(1,2)=T3*(XR*T6) BP(2,1)=T3*YR*T5 BP(2,2)=T3+T3*(T+YR*T6) AMP11=BP(1,1)*(PX*BP(1,1)+PXY*BP(1,2))+BP(1,2)*(PXY*BP(1,1)+PY* *BP(1,2)) AMP12=BP(1,1)*(PX*BP(2,1)+PXY*BP(2,2))+BP(1,2)*(PXY*BP(2,1)+PY* *BP(2,2)) AMP21=BP(2,1)*(PX*BP(1,1)+PXY*BP(1,2))+BP(2,2)*(PXY*BP(1,1)+PY* *BP(1,2)) AMP22=BP(2,1)*(PX*BP(2,1)+PXY*BP(2,2))+BP(2,2)*(PXY*BP(2,1)+PY*BP *(2,2)) DET=AMP11*AMP22-AMP12*AMP21 MP(1,1)=AMP22/DET MP(1,2)=-AMP12/DET MP(2,1)=-AMP21/DET MP(2,2)=AMP11/DET 600 W(1)=W1*MP(1,1)+W2*MP(1,2) W(2)=W1*MP(2,1)+W2*MP(2,2) DO 110 II=1,2 DO 120 K=1,3 MPA(II,K)=0.0 DO 120 JJ=1,2 120 MPA(II,K)=MPA(II,K)+MP(II,JJ)*AP2(JJ,K) 110 CONTINUE IF(IIO.EQ.0)GO TO 5150 1005 CONTINUE DO 1951 IK=1,3 KK=(KN-2)*3+IK DO 1952 IM=1,NUN UP1(IM)=0.0 DO 1952 IN=1,3 1952 UP1(IM)=UP1(IM)+NP21(IM,IN)*NP22(IN,KK) DO 1953 IL=1,NUN B(IL)=0.0 DO 1953 IP=1,NUN IOJ=((IP-1)**2+IP-1)/2+IL IF(IL.GT.IP)IOJ=((IL-1)**2+IL-1)/2+IP 1953 B(IL)=B(IL)+C(IOJ)*UP1(IP) DO 1957 KI=1,NUN 1957 NP21(KI,IK+3)=UP1(KI) AXX(IK)=0.0 DO 1954 IR=1,NUN 1954 AXX(IK)=AXX(IK)+B(IR)*UP1(IR) AXX(IK)=DSQRT(AXX(IK)) AXX(IK)=AXX(IK)*SEUW IF(IK.EQ.1)GO TO 1951 IF(IK.EQ.3)GO TO 1961 AXX(4)=0.0 DO 1959 ITJ=1,NUN 1959 AXX(4)=AXX(4)+B(ITJ)*NP21(ITJ,4) GO TO 1951 1961 AXX(5)=0.0 DO 1962 ITJ=1,NUN 1962 AXX(5)=AXX(5)+B(ITJ)*NP21(ITJ,4) AXX(6)=0.0 DO 1963 ITJ=1,NUN 1963 AXX(6)=AXX(6)+B(ITJ)*NP21(ITJ,5) 1951 CONTINUE CX=AXX(1) CY=AXX(2) CZ=AXX(3) SEGX=SEGX+CX SEGY=SEGY+CY SEGZ=SEGZ+CZ N=J-1 LL=(AL(N)+0.01)/1000000000 A1L=LL AR=AL(N)-A1L*1000000000.0+0.01 JSO=AR JSO=JSO/1000 AXX(4)=AXX(4)*(SEUW**2) AXX(5)=AXX(5)*(SEUW**2) AXX(6)=AXX(6)*(SEUW**2) WRITE(IWP,51) 51 FORMAT(1H1) WRITE(IWP,50) 50 FORMAT(//5X,'ST. DIV. OF ADJ. COORDS.'/5X,24('-')/16X,'PT.#',11X, *'X',14X,'Y',14X,'Z',8X,'CORR. XY',5X,'CORR. XZ',5X,'CORR. YZ'/) WRITE(IWP,1955)JSO,(AXX(KL),KL=1,6) 1955 FORMAT(I20,3F15.3,3D13.2) V(1,1)=AXX(1)**2 V(2,2)=AXX(2)**2 V(3,3)=AXX(3)**2 V(1,2)=AXX(4) V(1,3)=AXX(5) V(2,3)=AXX(6) V(2,1)=V(1,2) V(3,2)=V(2,3) V(3,1)=V(1,3) CALL EIG(V,CCL,QK,ST,WD) ST1=ST1+ST(1) ST2=ST2+ST(2) ST3=ST3+ST(3) IF(IZO.EQ.1)GO TO 1000 DO 363 I=1,JCS DO 363 K=1,NU 363 NP21(I,K)=0.0 5150 CONTINUE DO 150 II=1,NU I2=IX*NU-NU+II DO 160 K=1,3 DO 161 JJ=1,2 161 NP21(I2,K)=NP21(I2,K)+AP1(JJ,II)*MPA(JJ,K) 160 CONTINUE 150 CONTINUE 501 CONTINUE IF(J.EQ.NP)IZO=1 IF(J.EQ.NP)GO TO 2005 J=J+1 GO TO 100 2005 I1=IMIN(IX) KN=KN+1 GO TO 1005 1000 DF=JC-ICO DV=JC-IV CEGX=SEGX/DF CEGY=SEGY/DF CEGZ=SEGZ/DV ST1=ST1/DF ST2=ST2/DF ST3=ST3/DV WRITE(IWP,51) IIII=JC-ICO JJJJ=JC-IV WRITE(IWP,1001)CEGX,ST1,IIII,CEGY,ST2,IIII,CEGZ,ST3,JJJJ 1001 FORMAT(///5X,'MEANS OF ALL NON-CONTROL POINTS'/5X,31('*')//10X, *'STANDARD ERRORS (OBJECT)',5X,'SEMI-MAJOR AXES (E.E.)'/9X,26('-'), *3X,24('-')/5X,'X-',F19.4,15X,F10.4,I12/ * 5X,'Y-',F19.4,15X,F10.4,I12/ * 5X,'Z-',F19.4,15X,F10.4,I12/1H1) RETURN END SUBROUTINE UPEC(UP1,INU,NP21,UP2,NP22,Y,JC,NAM,JJS,VEC,L,IL2,INC, *NG,UG,IG,IGN,JGN) IMPLICIT REAL*8(A-H,O-Z) REAL*8 Y(JJS) REAL*8 UP1(NAM),UP2(JJS),NP22(3,JJS),NP21(JJS,17),VEC(JJS) *,NG(IGN),UG(JGN) COMMON NPHO(8),IMAX(8),IMIN(8) NU=17 JCS=INC NUN=INU NAN=NUN IF(IG.EQ.0)GO TO 222 II=3*IG DO 15 I=1,II 15 Y(I)=UG(I)+UP2(I) DO 20 JP=1,II VEC(JP)=0.0 DO 25 J=1,II JPOJ=((J-1)**2+J-1)/2+JP IF(JP.GT.J)JPOJ=((JP-1)**2+JP-1)/2+J 25 VEC(JP)=VEC(JP)+NG(JPOJ)*Y(J) 20 CONTINUE 222 IX=IG+1 DO 30 I=IX,JC DO 30 J=1,3 JP=3*I-3+J VEC(JP)=0.0 DO 31 K=1,3 KP=3*I-3+K 31 VEC(JP)=VEC(JP)+NP22(J,KP)*UP2(KP) 30 CONTINUE IF(L.EQ.4)RETURN NF=NUN/NU REWINDIL2 II=0 DO 40 J=1,NF IKX=IMAX(J) IKY=IMIN(J) DO 40 I=1,NU READ(IL2)(NP21(KK,I),KK=IKY,IKX) II=II+1 DO 50 K=IKY,IKX 50 UP1(II)=UP1(II)-NP21(K,I)*VEC(K) 40 CONTINUE RETURN END SUBROUTINE MORMAL(JC,JJS,NUM,NAM,NP22,NP11,PX1,NE,VEC,NP21,INC,INU *,IL2,C,MOM,IWP,ICAT,XP21,ITR,NITR,NPF,INOT,INPT,IG,NO,IGN,NU) IMPLICIT REAL*8(A-H,O-Z) DIMENSION ICAT(NAM),NPF(250),INOT(INPT) REAL*8 NP21(JJS,17),NP22(3,JJS),NP11(17,NUM),C(MOM),NE(NAM,17), *VEC(JJS),XP21(JJS,17),NO(IGN) REAL*4 PX1(8,17) COMMON NPHO(8),IMAX(8),IMIN(8) REWINDIL2 JCS=INC NUN=INU NAN=NUN IDR=0 IQ=0 IOP=0 DO 102 I=1,MOM 102 C(I)=0.0 IT=IG+1 NFI=NUN/NU L=1 IIJ=IMIN(1) JJJ=IMAX(1) DO 16 J=1,NU READ(IL2)(NP21(KK,J),KK=1,JJJ) 16 CONTINUE C C DO 10 K=1,NFI IDR=IDR+IQ IQ=NPF(K) JKL=(K-1)*NU DO 33 JH=1,NU JKLOJH=JKL+JH 33 ICAT(JKLOJH)=0 DO 11 I4=1,JCS DO 11 I5=1,NU 11 XP21(I4,I5)=0.0 IX=K KKX=IMAX(IX) IIX=IMIN(IX) I2=KKX/3 I1=(IIX-1)/3+1 333 IF(IG.EQ.0)GO TO 222 IF(I1.GT.IG)GO TO 222 DO 20 L=1,NU DO 20 J=1,IG DO 20 II=1,3 JJ=3*J-3+II XP21(JJ,L)=0.0 DO 15 I=I1,I2 IF(I.GT.IG)GO TO 20 I32=3*I-2 I31=3*I-1 I30=3*I JJJ=((I32-1)**2+I32-1)/2+JJ IF(JJ.GT.I32)JJJ=((JJ-1)**2+JJ-1)/2+I32 KKK=((I31-1)**2+I31-1)/2+JJ IF(JJ.GT.I31)KKK=((JJ-1)**2+JJ-1)/2+I31 LLL=((I30-1)**2+I30-1)/2+JJ IF(JJ.GT.I30)LLL=((JJ-1)**2+JJ-1)/2+I30 15 XP21(JJ,L)=XP21(JJ,L)+NP21(I32,L)*NO(JJJ)+NP21(I31,L)*NO(KKK)+ *NP21(I30,L)*NO(LLL) 20 CONTINUE 222 CONTINUE IF(I2.LE.IG)GO TO 351 DO 26 L=1,NU DO 25 J=I1,I2 IF(J.LT.IT)GO TO 25 DO 29 II=1,3 JJ=3*J-3+II 29 XP21(JJ,L)=NP21(3*J-2,L)*NP22(II,3*J-2)+NP21(3*J-1,L)* *NP22(II,3*J-1)+NP21(3*J,L)*NP22(II,3*J) 25 CONTINUE 26 CONTINUE 351 CONTINUE REWINDIL2 II=0 JDR=0 JQ=0 IL=K+1 DO 40 I=1,K JDR=JDR+JQ JQ=NPF(I) J1=IMIN(I) J2=IMAX(I) DO 40 J=1,NU II=II+1 READ(IL2)(NP21(M,J),M=J1,J2) DO 512 KG=1,NU NE(II,KG)=0.0 IF(III.GT.J2)GO TO 512 DO 91 JO=1,JQ ITIT=JDR+JO KO=INOT(ITIT) IF(KO.GT.I2.OR.KO.LT.I1)GO TO 91 DO 41 IS=1,3 KK=3*KO-3+IS 41 NE(II,KG)=NE(II,KG)-NP21(KK,J)*XP21(KK,KG) 91 CONTINUE 512 CONTINUE 40 CONTINUE IF(K.EQ.NFI)GO TO 401 J1=IMIN(K+1) J2=IMAX(IX+1) DO 82 J=1,NU 82 READ(IL2)(NP21(M,J),M=J1,J2) 401 CONTINUE DO 45 I=1,NU M=JKL+I DO 43 III=1,NU JJ=JKL+III 43 NE(JJ,I)=NP11(III,M)+NE(JJ,I) NE(M,I)=NE(M,I)+PX1(K,I) 45 CONTINUE DO 140 M=1,NU J=JKL+M IDX=0 IF(ITR.EQ.NITR)IDX=1 IFL=0 DO 120 I=1,J IF(NE(I,M).NE.0.0D0)GO TO 121 IFL=IFL+1 IF(IDX.EQ.0)ICAT(J)=IFL IF(IDX.EQ.0)GO TO 120 121 IDX=1 IOP=IOP+1 C(IOP)=NE(I,M) 120 CONTINUE IF(IOP.GT.MOM)WRITE(IWP,501) IF(IOP.GT.MOM)STOP 140 CONTINUE 10 CONTINUE 501 FORMAT(/'***ERROR*** TOO MANY PHOTOS'//) RETURN END SUBROUTINE QHLSKS(M,R,B,X,Y,ITR,N1,INU,Q,MOM,ICAT) C ********************************************************************* C *SUBROUTINE CHLSKS C * C * THIS SUBROUTINE SOLVES VECTOR EQUATION AX+B=0 FOR C * VECTOR X, WHERE VECTOR B IS GIVEN AND WHERE R, THE C * CHOLESKY DECOMPOSITION OF THE POSITIVE DEFINITE MATRIX C * A, IS KNOWN.. C * C ********************************************************************* C IMPLICIT REAL*8(A-H,O-Z) INTEGER Q,ICAT(M) REAL*8 X(M),Y(M) REAL*8 R(MOM),B(M) N=INU+N1 IF(ITR.LT.Q)N=N-N1 IPX=0 DO 30 K=1,N IPX=IPX+ICAT(K) S=B(K) IF(K-1)25,25,10 10 KLSS1=K-1 J=ICAT(K)+1 DO 20 I=J,KLSS1 IOK=((K-1)**2+K-1)/2+I-IPX S=S+R(IOK)*Y(I) 20 CONTINUE 25 KOK=((K-1)**2+K-1)/2+K-IPX Y(K)=-S/R(KOK) 30 CONTINUE DO 50 IDUM=1,N I=N+1-IDUM S=Y(I) IF(N-I)45,45,34 34 IPLS1=I+1 DO 40 K=IPLS1,N IKX=IKX+ICAT(K) IF(I.LE.ICAT(K))GO TO 40 IOK=((K-1)**2+K-1)/2+I-IKX S=S-R(IOK)*X(K) 40 CONTINUE 45 IOI=((I-1)**2+I-1)/2+I-IPX X(I)=S/R(IOI) IPX=IPX-ICAT(I) IKX=IPX 50 CONTINUE RETURN END SUBROUTINE CHECK(L,X,Y,Z,LC,XC,YC,ZC,IC,NB,INB,IN,AL,NCH,IWP) IMPLICIT REAL*8(A-H,O-Z) DIMENSION L(NB),LC(NCH) REAL*8 X(NB),Y(NB),Z(NB),XC(NCH),YC(NCH),ZC(NCH) REAL*8 AL(NB),AR,A1L NP=INB SX=0. SY=0.0 SZ=0.0 SEX=0.0 YEX=0.0 ZEX=0.0 WRITE(IWP,17) 17 FORMAT(//5X,'CHECK PT. DIFFS.'/5X,16('-')/) DO 30 I=1,IN J=LC(I) A=XC(I) B=YC(I) C=ZC(I) DO 980 II=1,NP LL=(AL(II)+0.01)/1000000000 A1L=LL AR=AL(II)-A1L*1000000000.0+0.01 JR=AR IR=JR/1000 980 IF(LC(I).EQ.IR)K=II AA=X(K)-A SEX=SEX+AA AD=AA**2 SX=SX+AD BB=Y(K)-B YEX=YEX+BB BD=BB**2 SY=SY+BD CC=Z(K)-C ZEX=ZEX+CC WRITE(IWP,22)J,AA,BB,CC CC=CC**2 SZ=SZ+CC 30 CONTINUE ZN=IN RMX=DSQRT(SX/ZN) RMY=DSQRT(SY/ZN) RMZ=DSQRT(SZ/ZN) SEX=SEX/ZN YEX=YEX/ZN ZEX=ZEX/ZN WRITE(IWP,21)RMX,RMY,RMZ WRITE(IWP,23)SEX,YEX,ZEX 23 FORMAT(//5X,'MEAN=',3F10.4) 22 FORMAT(I10,3F10.3) 21 FORMAT(///5X,'RMS =',3F10.4,2X,'( METERS )') RETURN END SUBROUTINE EIGENZ(T,R,VA,A,N,NX,NV) C IMPLICIT REAL *8(A-H,O-Z) DIMENSION T(NX,NX),VA(NX),A(6),R(9) C K=0 DO 10 J=1,N DO 10 I=1,J K=K+1 A(K)=T(I,J) 10 CONTINUE NN=N*(N+1)/2 20 IF (NV-1) 30,60,30 30 IQ=-N DO 50 J=1,N IQ=IQ+N DO 50 I=1,N IJ=IQ+I R(IJ)=0.0D0 IF (I-J) 50,40,50 40 R(IJ)=1.0D0 50 CONTINUE C C COMPUTE INITIAL AND FINAL NORMS ANORM AND ANORMX C 60 ANORM=0.0D0 DO 80 I=1,N DO 80 J=I,N IF (I-J) 70,80,70 70 IA=I+(J*J-J)/2 ANORM=ANORM+A(IA)*A(IA) 80 CONTINUE IF (ANORM) 370,370,90 90 ANORM=DSQRT(2.0D2*ANORM) ANRMX=ANORM*1.0D-06/DFLOAT(N) C C INITIALIZE INDICATORS AND COMPUTE THRESHOULD, THR C IND=0 THR=ANORM 100 THR=THR/DFLOAT(N) 110 L=1 120 M=L+1 C C COMPUTE SIN AND COS C 130 MQ=(M*M-M)/2 LQ=(L*L-L)/2 LM=L+MQ 140 IF (DABS(A(LM))-THR) 300,150,150 150 IND=1 LL=L+LQ MM=M+MQ X=0.5D0*(A(LL)-A(MM)) 160 Y=-A(LM)/DSQRT(A(LM)*A(LM)+X*X) IF (X) 170,180,180 170 Y=-Y 180 SINX=Y/DSQRT(2.0D0*(1.0D0+(DSQRT(1.0D0-Y*Y)))) SINX2=SINX*SINX 190 COSX=1.0D0-SINX2 COSX=DSQRT(COSX) COSX2=COSX*COSX SINCS=SINX*COSX C C ROTATE L AND M COLUMNS C ILQ=N*(L-1) IMQ=N*(M-1) DO 290 I=1,N IQ=(I*I-I)/2 IF (I-L) 200,270,200 200 IF (I-M) 210,270,220 210 IM=I+MQ GO TO 230 220 IM=M+IQ 230 IF (I-L) 240,250,250 240 IL=I+LQ GO TO 260 250 IL=L+IQ 260 X=A(IL)*COSX-A(IM)*SINX A(IM)=A(IL)*SINX+A(IM)*COSX A(IL)=X 270 IF (NV-1) 280,290,280 280 ILR=ILQ+I IMR=IMQ+I X=R(ILR)*COSX-R(IMR)*SINX R(IMR)=R(ILR)*SINX+R(IMR)*COSX R(ILR)=X 290 CONTINUE X=2.0D0*A(LM)*SINCS Y=A(LL)*COSX2+A(MM)*SINX2-X X=A(LL)*SINX2+A(MM)*COSX2+X A(LM)=(A(LL)-A(MM))*SINCS+A(LM)*(COSX2-SINX2) A(LL)=Y A(MM)=X C C TESTS FOR COMPLITION C C TEST FOR M LAST COLUMN C 300 IF (M-N) 310,320,310 310 M=M+1 GO TO 130 C C TEST FOR L SECOND FROM LAST COLUMN C 320 IF (L-(N-1)) 330,340,330 330 L=L+1 GO TO 120 340 IF (IND-1) 360,350,360 350 IND=0 GO TO 110 C C COMPARE THRESHOULD WITH FINAL NORM C 360 IF (THR-ANRMX) 370,370,100 C C SORT EIGENVALUES AND EIGENVECTORS C 370 IQ=-N DO 410 I=1,N IQ=IQ+N LL=I+(I*I-I)/2 JQ=N*(I-2) DO 410 J=I,N JQ=JQ+N MM=J+(J*J-J)/2 IF (A(LL)-A(MM)) 380,410,410 380 X=A(LL) A(LL)=A(MM) A(MM)=X IF (NV-1) 390,410,390 390 DO 400 K=1,N ILR=IQ+K IMR=JQ+K X=R(ILR) R(ILR)=R(IMR) 400 R(IMR)=X 410 CONTINUE K=0 DO 420 I=1,N K=K+I 420 VA(I)=A(K) RETURN END SUBROUTINE EIG(A,CL,CK,S,WD) IMPLICIT REAL*8(A-H,O-Z) C C PROGRAM FOR THE ANALYSIS OF ERROR ELLIPSOID C DIMENSION A(3,3),R(9),EV(3),WS(6),S(3),B(3,3) INTEGER*4 XYZ(3)/'X','Y','Z'/ EQUIVALENCE (R(1),B(1)) C C READ COVARIANCE MATRIX C C SUMA=0.0 DO 50 I=1,3 50 SUMA=SUMA+A(I,I) C WRITE(6,11) WRITE(6,12) ((A(I,J),J=1,3),I=1,3) C C EIGENVALUES AND EIGENVECTORS C CALL EIGENZ(A,R,EV,WS,3,3,0) C SUMS=0.0 NEG=1 DO 60 I=1,3 IF(EV(I).LT.0.0) NEG=-1 SUMS=SUMS+EV(I) MII=4-I S(I)=DSQRT(DABS(EV(MII)))*CK DO 60 J=1,3 JKLL=9-I*3+J 60 A(I,J)=R(JKLL) C CCCC WRITE(6,13) (EV(4-I),(A(I,J),J=1,3),I=1,3) C C DO 65 I=1,3 C DO 65 J=1,3 C B(I,J)=0.0 C DO 65 K=1,3 C 65 B(I,J)=B(I,J)+A(K,I)*A(K,J) C PAI=3.1415926535898 DO 70 I=1,3 DO 70 J=1,3 70 A(I,J)=DARCOS(A(I,J))*180.0D0/PAI C IF(NEG.EQ.-1) GO TO 100 C CCCC WRITE(6,14)WD,CL CCCC WRITE(6,17) CCCC WRITE(6,15) (XYZ(I),S(I),XYZ(I),(A(I,J),J=1,3),I=1,3) C WRITE(6,16) ((B(I,J),J=1,3),I=1,3) GO TO 110 C CC100 WRITE(6,18) 110 CONTINUE C 10 FORMAT(6G13.3) 11 FORMAT(10X,1X/,T35,'*********************************',/, * T35,'* *',/, * T35,'* ANALYSIS OF ERROR-ELLIPSOID *',/, * T35,'* *',/, * T35,'*********************************',///) 12 FORMAT(1H0, 5X, '(1) COVARIANCE MATRIX(OBJECT):',//,(T10,3D16.6)) 13 FORMAT(1H0,/,6X,'(2) EIGENVALUES AND EIGENVECTORS OF COVARIANCE ', * 'MATRIX:',// ,T14,'EIGENVALUES',T41,'EIGENVECTORS (DIRECTI', * 'ON COSINES)',//,(T9, F15.6,T34,F15.6,2F12.6)) 14 FORMAT(1H0,/, 6X,'(3) PARAMETERS OF ',A8,' ERROR-ELLIPSOID (', * F4.1,'% CONFIDENCE REGION):') 15 FORMAT(1H0,T41,'SPATIAL ANGLES (IN DEGREES) BETWEEN AXES',//, * T14,'SEMI-MAJOR AXES', T53,'X(OBJECT)',2X,'Y(OBJECT)',2X, * 'Z(OBJECT)',//,(T14,A1,' =',F12.6,T41,A1,'(ELLIPSOID)', * F7.2,2F10.2)) C 16 FORMAT(1H0,/, 6X,'(4) ORTHOGONALITY TEST',11H (M'M = I):,//, C * (T10,3F16.6)) 17 FORMAT(1H+,T25,'________') 18 FORMAT(1H0,//,6X,'***(ERROR)*** NEGATIVE EIGENVALUE(S) DETECTED', * //,20X,'CHECK INPUT DATA FOR ROUND-OFF ERRORS') C RETURN END SUBROUTINE COEFS(C,A,B,QVX,QVY,MOM,IX,X,JJS,N,A2,KN,NP21) IMPLICIT REAL*8(A-H,O-Z) REAL*8 C(MOM),A(2,17),B(2,2),X(JJS),Q(2),N(3,JJS),A2(2,3) *,NP21(JJS,17),M12(17) NU=17 IT=(IX-1)*NU IS=(KN-1)*3 DO 20 I=1,2 DO 10 K=1,NU KK=IT+K X(K)=0.0 DO 10 J=1,NU ITOJ=IT+J KOJ=((ITOJ-1)**2+ITOJ-1)/2+KK IF(KK.GT.ITOJ)KOJ=((KK-1)**2+KK-1)/2+ITOJ 10 X(K)=X(K)+C(KOJ)*A(I,J) AX=0.0 DO 15 II=1,NU 15 AX=AX+A(I,II)*X(II) AOR=AX*B(I,I) DO 40 K=1,3 KK=IS+K X(K)=0.0 DO 40 J=1,3 KJ=IS+J 40 X(K)=X(K)+N(K,KJ)*A2(I,J) DO 45 II=1,3 45 AX=AX+A2(I,II)*X(II) ACO=AX*B(I,I) ACO=ACO-AOR DO 50 II=1,NU M12(II)=0.0 DO 50 JJ=1,3 JK=IS+JJ 50 M12(II)=M12(II)+NP21(JK,II)*X(JJ) DO 60 II=1,NU I1=IT+II X(II)=0.0 DO 60 JJ=1,NU J1=IT+JJ IOJ=((J1-1)**2+J1-1)/2+I1 IF(I1.GT.J1)IOJ=((I1-1)**2+I1-1)/2+J1 60 X(II)=X(II)+C(IOJ)*M12(JJ) DO 65 II=1,NU 65 AX=AX+M12(II)*X(II) ADO=AX*B(I,I) ADO=ADO-ACO-AOR DO 70 II=1,NU 70 AX=AX-2.0*A(I,II)*X(II) 20 Q(I)=1.0-AX*B(I,I) QVX=Q(1) QVY=Q(2) RETURN END SUBROUTINE SYINV3(A) IMPLICIT REAL*8(A-H,O-Z) REAL*8 A(3,3) DET=A(1,1)*A(2,2)-A(1,2)**2 PX=A(2,2)/DET PY=A(1,1)/DET PXY=-A(1,2)/DET A1=A(1,3)*PX+A(2,3)*PXY A2=A(1,3)*PXY+A(2,3)*PY EP=A(3,3)-A1*A(3,1)-A2*A(3,2) A(1,1)=PX+A1/EP*A1 A(1,2)=PXY+A1/EP*A2 A(2,1)=PXY+A2/EP*A1 A(2,2)=PY+A2/EP*A2 A(1,3)=-A1/EP A(2,3)=-A2/EP A(3,1)=A(1,3) A(3,2)=A(2,3) A(3,3)=1.0/EP RETURN END SUBROUTINE SMINV(A,N1,NN,SE,MOM) IMPLICIT REAL*8(A-H,O-Z) REAL*8 SE(N1),A(MOM) DET=A(1)*A(3)-A(2)**2 PX=A(3)/DET PY=A(1)/DET PXY=-A(2)/DET A1=A(4)*PX+A(5)*PXY A2=A(4)*PXY+A(5)*PY EP=A(6)-A1*A(4)-A2*A(5) A(1)=PX+A1/EP*A1 A(2)=PXY+A1/EP*A2 A(3)=PY+A2/EP*A2 A(4)=-A1/EP A(5)=-A2/EP A(6)=1.0/EP N=3 200 N=N+1 M=N-1 DO 11 I=1,M SE(I)=0.0D00 DO 10 J=1,M JON=((N-1)**2+N-1)/2+J IOJ=((J-1)**2+J-1)/2+I IF(I.GT.J)IOJ=((I-1)**2+I-1)/2+J 10 SE(I)=SE(I)+A(IOJ)*A(JON) 11 CONTINUE EP=0.0D00 DO 12 J=1,M NOJ=((N-1)**2+N-1)/2+J 12 EP=EP+A(NOJ)*SE(J) NON=((N-1)**2+N-1)/2+N EB=A(NON)-EP EI=1.0D00/EB DO 100 I=1,M DO 100 J=I,M IOJ=((J-1)**2+J-1)/2+I 100 A(IOJ)=A(IOJ)+SE(I)*SE(J)*EI DO 150 I=1,M ION=((N-1)**2+N-1)/2+I 150 A(ION)=-SE(I)*EI A(NON)=EI IF(NN.GT.N)GO TO 200 RETURN END SUBROUTINE UTINV(A,N1,NN,SE,MOM) IMPLICIT REAL*8(A-H,O-Z) REAL*8 A(MOM),SE(N1) A(1)=1.0/A(1) A(2)=-A(2)*A(1)/A(3) A(3)=1.0/A(3) N=2 200 M=N N=N+1 DO 10 I=1,M SE(I)=0.0 DO 11 J=I,M JON=((N-1)**2+N-1)/2+J IOJ=((J-1)**2+J-1)/2+I 11 SE(I)=SE(I)+A(IOJ)*A(JON) 10 CONTINUE NON=((N-1)**2+N-1)/2+N DO 12 I=1,M ION=((N-1)**2+N-1)/2+I 12 A(ION)=-SE(I)/A(NON) A(NON)=1.0/A(NON) IF(NN.GT.N)GO TO 200 RETURN END SUBROUTINE QHLSK(M,A,R,INDEF,ITR,N1,INU,Q,MOM,IWP,ICAT) C C ********************************************************************* C *SUBROUTINE CHLSK C * C * THIS SUBROUTINE COMPUTES THE UPPER TRIANGULAR MATRIX R C * WHICH, WHEN MULTIPLIED BY ITS OWN TRANSPOSE, YIELDS C * THE GIVEN SYMMETRIC MATRIX A, BY THE CHOLESKY METHOD C * C * C ********************************************************************* C IMPLICIT REAL*8(A-H,O-Z) INTEGER Q,ICAT(M) REAL*8 A(MOM),R(1,M) N=INU+N1 IF(ITR.LT.Q)N=N-N1 INDEF=0 DO 30 IP=1,N INDEF=INDEF+ICAT(IP) IPOIP=((IP-1)**2+IP-1)/2+IP-INDEF IF(A(IPOIP))05,05,10 5 WRITE(IWP,6)IP 6 FORMAT(/5X,'***ERROR*** ARRAY NOT POSITIVE DEFINITE AT COMP#',I4) STOP 10 R(01,IP)=DSQRT(A(IPOIP)) IF(IP-N)12,40,40 12 IPPLS1=IP+1 IPX=INDEF DO 20 K=IPPLS1,N IPX=IPX+ICAT(K) IF(IP.LE.ICAT(K))R(1,K)=0.0D0 IF(IP.LE.ICAT(K))GO TO 20 IPOK=((K-1)**2+K-1)/2+IP-IPX R(01,K)=A(IPOK)/R(01,IP) 20 CONTINUE IGX=INDEF DO 31 I=IPPLS1,N IPX=IGX IGX=IGX+ICAT(I) DO 31 K=I,N IPX=IPX+ICAT(K) IF(I.LE.ICAT(K))GO TO 31 IOK=((K-1)**2+K-1)/2+I-IPX A(IOK)=A(IOK)-R(01,I)*R(01,K) 31 CONTINUE IPX=INDEF-ICAT(IP) IF(IP.EQ.1)IPX=0 DO 50 IX=IP,N IPX=IPX+ICAT(IX) IF(IP.LE.ICAT(IX))GO TO 50 IPIX=((IX-1)**2+IX-1)/2+IP-IPX A(IPIX)=R(1,IX) 50 CONTINUE 30 CONTINUE 40 DO 60 IX=IP,N IPX=INDEF IPIX=((IX-1)**2+IX-1)/2+IP-IPX A(IPIX)=R(1,IX) 60 CONTINUE RETURN END SUBROUTINE KHLSK(M,A,R,MOM,N,IWP) IMPLICIT REAL*8(A-H,O-Z) REAL*8 A(MOM),R(1,M) DO 30 IP=1,N IPOIP=((IP-1)**2+IP-1)/2+IP IF(A(IPOIP))05,05,10 5 WRITE(IWP,6)IP 6 FORMAT(/5X,'***ERROR*** ARRAY NOT POSITIVE DEFINITE AT COMP#',I4, *'(GEOD. NOR. EQ.)') STOP 10 R(01,IP)=DSQRT(A(IPOIP)) IF(IP-N)12,40,40 12 IPPLS1=IP+1 DO 20 K=IPPLS1,N IPOK=((K-1)**2+K-1)/2+IP R(01,K)=A(IPOK)/R(01,IP) 20 CONTINUE DO 31 I=IPPLS1,N DO 31 K=I,N IOK=((K-1)**2+K-1)/2+I A(IOK)=A(IOK)-R(01,I)*R(01,K) 31 CONTINUE DO 50 IX=IP,N IPIX=((IX-1)**2+IX-1)/2+IP A(IPIX)=R(1,IX) 50 CONTINUE 30 CONTINUE 40 DO 60 IX=IP,N IPIX=((IX-1)**2+IX-1)/2+IP A(IPIX)=R(1,IX) 60 CONTINUE RETURN END SUBROUTINE DNORM(NG,UG,NAN,NP,XG,YG,ZG,IG,ITR,S,WT,INB,IGM,L,XO, *IGN,AG,WG,AL,JGN,IRC,IWP,ICOD) IMPLICIT REAL*8(A-H,O-Z) REAL*4 S(IGM,IGM),WT(IGM,IGM) REAL*8 UG(JGN),XG(NP),YG(NP),ZG(NP),NG(IGN),XO(JGN) REAL*8 AL(NP),AR,A1L,AG(1,JGN) DIMENSION L(NP) 10 FORMAT(2I10,F10.3,F10.3) IF(ITR.EQ.1)WRITE(IWP,11) 11 FORMAT(//5X,'READ IN HEIGHT DIFFERENCES'/11X,'FROM',8X,'TO',10X, *'ST.E.',10X,'HEIGHT DIFF.'/) 111 FORMAT(//5X,'HT. DIFF. FROM ADJ. COORDS'/11X,'FROM',8X,'TO',10X, *'DIFF.',10X,'HEIGHT DIFF.'/) 12 FORMAT(5X,2I10,5X,F10.3,5X,F10.3) IGS=IG*3 IF(ICOD.EQ.12)GO TO 333 DO 21 II=1,IGS 21 UG(II)=0.0 DO 20 JJ=1,IGN 20 NG(JJ)=0.0 333 CONTINUE 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)NOD DO 77 KM=1,NOD READ(IRC,10)I,J,P,SS WRITE(IWP,12)I,J,P,SS 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 77 CONTINUE WRITE(IWP,13)IG,NOD 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,NOD 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=XO(J3)-XO(I3) 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 AG(1,I3-2)=0.0 AG(1,J3-2)=0.0 AG(1,I3)=-1.0*WT(I,J) AG(1,J3)=1.0*WT(I,J) AG(1,J3-1)=0.0 AG(1,I3-1)=0.0 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 IF(I3OIA.GT.I3OJA)GO TO 301 INDL=((I3OJA-1)**2+I3OJA-1)/2+I3OIA NG(INDL)=NG(INDL)+AG(1,I3OIA)*AG(1,I3OJA)/WT(I,J) 301 CONTINUE DO 302 IB=1,3 IA=IB-1 DO 302 JB=1,3 JA=JB-1 J3OJA=J3-JA I3OIA=I3-IA INDL=((J3OJA-1)**2+J3OJA-1)/2+I3OIA IF(I3OIA.GT.J3OJA)INDL=((I3OIA-1)**2+I3OIA-1)/2+J3OJA NG(INDL)=NG(INDL)+AG(1,I3OIA)*AG(1,J3OJA)/WT(I,J) 302 CONTINUE DO 304 IB=1,3 IA=IB-1 DO 304 JB=1,3 JA=JB-1 J3OJA=J3-JA J3OIA=J3-IA IF(J3OIA.GT.J3OJA)GO TO 304 INDL=((J3OJA-1)**2+J3OJA-1)/2+J3OIA NG(INDL)=NG(INDL)+AG(1,J3OIA)*AG(1,J3OJA)/WT(I,J) 304 CONTINUE 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 SUBROUTINE NORM2(NP22,NG,PX2,JJS,JC,IG,SE,M,NAM,IRM,INC,INU,IUY, *IWP,IL1,R) IMPLICIT REAL*8(A-H,O-Z) REAL*4 PX2(1,JJS) REAL*8 NG(IUY),R(1,NAM) REAL*8 NP22(3,JJS),SE(JJS),M(3,3) NM=3*IG JCS=INC NUN=INU DO 100 I=1,JC III=3*I NP22(1,III-2)=NP22(1,III-2)+PX2(1,III-2) NP22(2,III-1)=NP22(2,III-1)+PX2(1,III-1) 100 NP22(3,III-0)=NP22(3,III-0)+PX2(1,III) IF(IG.EQ.0)GO TO 1000 DO 200 I=1,IG III=3*I KII=III+2 JK=((KII-5)**2+KII-5)/2+KII-4 NG(JK)=NG(JK)+NP22(1,III-2) JK=((KII-4)**2+KII-4)/2+KII-4 NG(JK)=NG(JK)+NP22(1,III-1) JK=((KII-3)**2+KII-3)/2+KII-4 NG(JK)=NG(JK)+NP22(1,III) JK=((KII-4)**2+KII-4)/2+KII-3 NG(JK)=NG(JK)+NP22(2,III-1) JK=((KII-3)**2+KII-3)/2+KII-3 NG(JK)=NG(JK)+NP22(2,III) JK=((KII-3)**2+KII-3)/2+KII-2 NG(JK)=NG(JK)+NP22(3,III) 200 CONTINUE IF(IG.LE.20)GO TO 150 CALL KHLSK(NAM,NG,R,IUY,NM,IWP) CALL UTINV(NG,JJS,NM,SE,IUY) REWINDIL1 DO 10 I=1,NM DO 11 K=I,NM SE(K)=0.0 DO 11 J=K,NM IOJ=((J-1)**2+J-1)/2+I KOJ=((J-1)**2+J-1)/2+K 11 SE(K)=SE(K)+NG(IOJ)*NG(KOJ) WRITE(IL1)(SE(II),II=I,NM) 10 CONTINUE ENDFILEIL1 REWINDIL1 DO 20 I=1,NM READ(IL1)(SE(II),II=I,NM) DO 21 J=I,NM IOJ=((J-1)**2+J-1)/2+I 21 NG(IOJ)=SE(J) 20 CONTINUE GO TO 1000 150 CALL SMINV(NG,JJS,NM,SE,IUY) 1000 CONTINUE K=1 DO 300 I=K,JC III=3*I DO 310 J=1,3 M(J,1)=NP22(J,III-2) M(J,2)=NP22(J,III-1) 310 M(J,3)=NP22(J,III) CALL SYINV3(M) DO 320 L=1,3 NP22(L,III-2)=M(L,1) NP22(L,III-1)=M(L,2) 320 NP22(L,III)=M(L,3) 300 CONTINUE RETURN END SUBROUTINE SNORM(NG,UG,NAN,NP,XG,YG,ZG,IG,ITR,S,WT,INB,IGM,L,XO, *IGN,AG,WG,AL,JGN,IRC,IWP) IMPLICIT REAL*8(A-H,O-Z) REAL*4 S(IGM,IGM),WT(IGM,IGM) REAL*8 UG(JGN),XG(NP),YG(NP),ZG(NP),NG(IGN),XO(JGN) REAL*8 AL(NP),AR,A1L,AG(1,JGN) DIMENSION L(NP) 10 FORMAT(2I10,F10.3,F10.3) 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 21 II=1,IGS 21 UG(II)=0.0 DO 20 JJ=1,IGN 20 NG(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 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 C NOO=4 C IF(KK.GT.NOO) GO TO 702 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 C 702 II=I C JJ=J C IF(ITR.EQ.0)GO TO 100 C AG(1,I3-2)=100.D0 C AG(1,I3-1)=100.D0 C AG(1,I3)=100.D0 C AG(1,J3-2)=0.D0 C AG(1,J3-1)=0.D0 C AG(1,J3)=0.D0 C 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 IF(I3OIA.GT.I3OJA)GO TO 301 INDL=((I3OJA-1)**2+I3OJA-1)/2+I3OIA NG(INDL)=NG(INDL)+AG(1,I3OIA)*AG(1,I3OJA)/WT(I,J) 301 CONTINUE DO 302 IB=1,3 IA=IB-1 DO 302 JB=1,3 JA=JB-1 J3OJA=J3-JA I3OIA=I3-IA INDL=((J3OJA-1)**2+J3OJA-1)/2+I3OIA IF(I3OIA.GT.J3OJA)INDL=((I3OIA-1)**2+I3OIA-1)/2+J3OJA NG(INDL)=NG(INDL)+AG(1,I3OIA)*AG(1,J3OJA)/WT(I,J) 302 CONTINUE DO 304 IB=1,3 IA=IB-1 DO 304 JB=1,3 JA=JB-1 J3OJA=J3-JA J3OIA=J3-IA IF(J3OIA.GT.J3OJA)GO TO 304 INDL=((J3OJA-1)**2+J3OJA-1)/2+J3OIA NG(INDL)=NG(INDL)+AG(1,J3OIA)*AG(1,J3OJA)/WT(I,J) 304 CONTINUE 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 SUBROUTINE DMSRAD(IDEG,IMIN,SEC,RAD) IMPLICIT REAL *8(A-H,O-Z) DEG=IABS(IDEG)+IABS(IMIN)/60.D0+DABS(SEC)/3600.D0 RAD=DEG*3.141592653589793D0/180.D0 RAD=RAD*ISIGN(1,IDEG) IF(IDEG.EQ.0)RAD=RAD*ISIGN(1,IMIN) IF(IDEG.EQ.0.AND.IMIN.EQ.0)RAD=RAD*DSIGN(1.D0,SEC) RETURN END