//GEBATV JOB ' , ',' ',MSGLEVEL=(1,1) /*JOBPARM S=119,R=2048,L=25 // EXEC FORTVCLG,PARM.FORT='LANGLVL(66)',RC=2048K,RL=2048K,RG=2048K //FORT.SYSIN DD * 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.....................4 * C * -MAX. NO. OF IMAGE PTS................320 * C * -MAX. NO. OF OBJECT PTS................80 * C * -MAX. NO. OF GEODETIC PTS..............15 * C * -UNLIMITED ; PTS./PHOTO ,CONTROL PTS. * C * & NO. OF DISTANCES&HEIGHTS * C * * C ************************************************************* C C IMPLICIT REAL*8(A-H,O-Z) DIMENSION L(321),LG(321),NF(080),NV(080),LC(080),LQ(080), *ICAT(68),INOT(321) REAL*4 XP(321),YP(321),PX1(17),PX2(1,240) *,X1(15,15),X2(15,15),X3(15,15),X4(15,15) REAL*8 XG(321),YG(321),ZG(321),XC(080),YC(080),ZC(080),XQ(080), *YQ(080),ZQ(080),X(68),B(68),P(240) *,UG(045),XO(045),NOG(1035),AG(1,045) REAL*8 AP1(2,17),MP(2,2),MPA(2,17),W(2),UP1(68),XP21(240,17), *NP11(17,68),M(3,3),BP(2,2),AP2(2,3),NP21(240,17),NP22(3,240), *UP2(240),NE(68,17),VEC(240),AL(321),R(1,68),C(2346) 3,FF(10),PF(10) COMMON NPHO(4),IMAX(4),IMIN(4) DATA NOB,NCO,NCH,IL1/80,80,80,21/ DATA IRC,IWP,IL2,IUY,INPT,IUNK,ITUN,IPTO/5,6,23,2346,321,68,68, *240/ DATA IGM,IGN,JGN/15,1035,45/ 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(4),INOT(INPT) REAL*4 XP(INPT),YP(INPT),PX1(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(4),IMAX(4),IMIN(4) 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 * 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 READ(IRC,201)(PX1(I),I=12,17),(PX1(J),J=1,3),P8 DO 1106 IK=4,11 1106 PX1(IK)=100.0 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 7414 PX1(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(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(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(4),IMAX(4),IMIN(4) 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(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(4),IMAX(4),IMIN(4) 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(4),IMAX(4),IMIN(4) 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(17) COMMON NPHO(4),IMAX(4),IMIN(4) 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(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 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 WRITE(6,14)WD,CL WRITE(6,17) 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 100 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 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 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 //GO.FT05F001 DD * 0 0 00 12 6700.0 7 3 1 0 (I10/2F10.3,2I10) (10X,I9,F16.3,F15.3) 5078 -297.886 300.093 1 1 11 -329.275 269.096 21 -329.185 263.979 31 -327.729 260.317 41 -325.161 255.875 51 -325.912 251.333 61 -318.639 251.603 71 -325.159 246.027 81 -323.035 240.682 91 -323.281 273.575 101 -305.241 270.444 111 -301.583 267.536 121 -300.628 265.497 131 -299.803 263.100 141 -297.396 260.071 151 -294.951 253.272 161 -295.219 245.998 171 -293.017 241.495 181 -279.118 276.397 191 -275.549 265.835 201 -272.161 255.449 211 -269.329 250.652 221 -268.415 239.722 231 -301.612 270.673 241 -290.169 266.054 251 -288.088 262.411 261 -284.362 257.993 271 -284.839 244.531 281 -274.589 290.463 291 -259.715 273.949 301 -248.963 299.009 311 -245.833 286.801 321 -333.767 275.930 331 -331.340 284.412 341 -310.008 278.622 351 -312.096 284.847 361 -305.931 286.055 371 -294.989 285.494 381 -289.495 283.788 391 -284.537 282.062 401 -298.159 298.260 411 -302.913 319.039 421 -257.858 321.037 431 -260.934 316.636 441 -258.173 308.340 1001 -261.511 341.061 1010 -292.181 297.660 1020 -308.588 272.585 1030 -330.067 272.230 1040 -284.771 277.016 11111 -219.729 198.264 70011 -325.567 340.437 70021 -338.485 325.571 70031 -337.742 303.772 70041 -334.915 290.266 70051 -308.770 298.006 70061 -276.559 291.396 70071 -250.887 309.493 70081 -210.873 362.914 70090 -262.155 340.871 70100 -251.697 347.751 70111 -364.740 221.157 70121 -281.179 230.774 70131 -293.204 313.567 90400 -241.519 199.727 666666 5079 -296.630 300.356 11 -381.594 269.307 21 -380.138 264.174 31 -377.529 260.491 41 -373.488 256.034 51 -373.209 251.519 61 -365.066 251.739 71 -370.921 246.207 81 -367.050 240.851 91 -376.542 273.723 101 -357.049 270.422 111 -351.811 267.511 121 -350.214 265.441 131 -348.717 263.033 141 -345.228 260.010 151 -341.214 253.172 161 -340.437 245.901 171 -336.845 241.377 181 -329.404 276.184 191 -322.476 265.577 201 -316.565 255.173 211 -312.343 250.340 221 -308.856 239.420 231 -352.567 270.654 241 -338.394 265.928 251 -335.167 262.263 261 -329.855 257.795 271 -327.844 244.363 281 -325.377 290.292 291 -305.410 273.585 301 -296.064 298.614 311 -290.358 286.387 321 -386.487 276.218 331 -383.872 284.687 341 -363.868 278.657 351 -366.620 284.917 361 -360.993 286.104 371 -349.308 285.400 381 -342.675 283.667 391 -336.806 281.885 401 -354.563 298.243 411 -359.915 319.169 421 -310.718 320.755 431 -313.764 316.380 441 -308.946 308.011 1001 -319.619 340.897 1010 -348.221 297.579 1020 -361.773 272.589 1030 -383.381 272.444 1040 -336.027 276.839 11111 -253.457 197.527 70011 -376.647 340.899 70021 -390.051 326.085 70031 -390.153 304.185 70041 -387.480 290.591 70051 -364.495 298.089 70061 -327.992 291.192 70071 -300.015 309.137 70081 -270.974 362.436 70090 -320.204 340.717 70100 -309.013 347.546 70111 -406.035 221.775 70121 -321.005 230.584 70131 -351.259 313.556 90400 -275.435 199.210 666666 5130 -298.737 300.675 11 -326.910 273.606 21 -327.110 268.612 31 -325.897 264.936 41 -323.585 260.441 51 -324.599 256.089 61 -317.391 255.805 71 -324.162 250.902 81 -322.413 245.610 91 -320.694 277.503 101 -302.853 272.995 111 -299.297 269.918 121 -298.511 267.835 131 -297.831 265.405 141 -295.621 262.291 151 -293.654 255.421 161 -294.456 248.320 171 -292.548 243.750 181 -275.944 276.917 191 -273.111 266.210 201 -270.469 255.694 211 -267.921 250.757 221 -267.835 239.990 231 -299.183 272.985 241 -287.907 267.596 251 -286.077 263.843 261 -282.613 259.219 271 -284.100 246.069 281 -270.197 290.737 291 -256.237 273.142 301 -242.927 297.661 311 -240.738 285.097 321 -330.857 280.699 331 -327.914 288.979 341 -307.137 281.497 351 -308.752 287.778 361 -302.504 288.530 371 -291.455 287.112 381 -285.991 285.035 391 -281.078 282.937 401 -293.761 300.135 411 -297.028 321.577 421 -250.400 320.659 431 -253.945 316.339 441 -251.825 307.705 1001 -252.559 341.421 1010 -287.727 299.092 1020 -306.093 275.331 1030 -327.515 276.717 1040 -281.688 277.944 11111 -221.675 194.841 70011 -318.205 345.232 70021 -332.204 330.916 70031 -332.955 308.841 70041 -331.070 295.093 70051 -304.489 300.720 70061 -272.168 291.787 70071 -244.065 308.429 70081 -196.604 361.179 70090 -253.236 341.273 70100 -241.633 347.821 70111 -364.023 230.617 70121 -281.404 232.387 70131 -287.589 315.241 90400 -243.794 198.504 666666 5131 -293.731 303.049 11 -372.431 272.345 21 -371.366 267.364 31 -369.110 263.723 41 -365.469 259.246 51 -365.569 254.941 61 -357.594 254.612 71 -363.770 249.768 81 -360.413 244.529 91 -367.071 276.122 101 -348.039 271.460 111 -343.157 268.414 121 -341.721 266.328 131 -340.459 263.911 141 -337.266 260.833 151 -333.869 253.971 161 -333.691 246.935 171 -330.527 242.403 181 -320.108 275.158 191 -314.141 264.581 201 -309.191 254.175 211 -305.452 249.270 221 -303.013 238.596 231 -343.607 271.457 241 -329.988 266.044 251 -327.124 262.346 261 -322.198 257.897 271 -321.398 244.715 281 -314.885 288.896 291 -296.430 271.379 301 -284.786 295.607 311 -280.223 283.213 321 -376.776 279.493 331 -373.623 287.775 341 -354.120 279.887 351 -356.392 286.223 361 -350.709 286.857 371 -339.123 285.351 381 -332.686 283.250 391 -327.001 281.129 401 -343.300 298.315 411 -347.133 319.784 421 -297.501 318.308 431 -300.960 314.081 441 -296.899 305.563 1001 -304.581 338.728 1010 -337.015 297.193 1020 -352.515 273.761 1030 -373.913 275.407 1040 -326.619 276.204 11111 -251.203 193.380 70011 -362.687 344.072 70021 -377.015 329.983 70031 -378.522 307.810 70041 -376.771 293.971 70051 -353.283 299.083 70061 -317.447 289.933 70071 -287.811 306.267 70081 -252.270 357.167 70090 -305.206 338.596 70100 -293.239 344.930 70111 -399.781 229.871 70121 -315.851 231.093 70131 -338.771 313.237 90400 -273.203 197.152 -66 (I6,F11.3,F14.3,F13.3,3F9.6) 11 1538.500 1429.412 1872.197 21 1506.290 1425.796 1849.594 31 1483.081 1415.958 1829.925 41 1454.168 1398.898 1803.592 51 1421.256 1399.133 1783.715 61 1427.092 1358.649 1768.094 71 1380.669 1390.318 1753.016 81 1336.917 1372.938 1714.973 91 1571.521 1401.983 1887.798 101 1572.544 1306.301 1867.637 111 1556.988 1282.709 1842.918 121 1544.684 1274.603 1831.380 131 1530.024 1266.439 1819.711 141 1511.734 1247.939 1799.880 151 1467.858 1222.357 1768.947 161 1417.795 1211.910 1745.949 171 1383.283 1188.679 1715.717 181 1632.468 1173.245 1848.821 191 1569.788 1128.816 1787.509 201 1502.561 1083.307 1735.051 211 1469.534 1051.300 1704.604 221 1381.654 1011.739 1640.227 231 1575.761 1287.218 1854.630 241 1556.389 1215.675 1809.472 251 1534.116 1196.649 1787.903 261 1506.376 1164.460 1755.953 271 1408.103 1139.390 1699.590 281 1713.595 1167.853 1860.243 291 1635.175 1042.248 1768.631 301 1794.400 1022.874 1803.888 311 1729.534 970.568 1751.503 321 1569.982 1459.273 1879.562 331 1616.101 1457.155 1878.323 341 1613.972 1340.060 1898.851 351 1643.862 1356.914 1909.674 361 1657.642 1327.654 1918.167 371 1666.302 1272.699 1909.311 381 1662.853 1241.997 1893.511 391 1658.679 1213.175 1880.364 401 1726.114 1303.324 1939.151 411 1819.169 1348.218 1949.967 421 1887.917 1124.845 1900.701 431 1862.329 1135.059 1899.178 441 1827.448 1103.390 1866.354 1001 1968.350 1172.800 1973.650 1010 1729.160 1274.240 1935.650 1020 1583.440 1325.790 1888.610 1030 1556.910 1435.230 1886.820 1040 1631.070 1204.760 1863.880 11111 1067.340 482.851 1426.256 70011 1914.790 1496.080 1868.398 70021 1820.912 1544.787 1871.013 70031 1708.358 1512.688 1879.470 70041 1642.146 1482.256 1879.309 70051 1713.218 1354.336 1928.841 70061 1716.139 1181.109 1870.426 70071 1846.042 1058.466 1841.405 70081 2066.540 1071.700 1985.350 70090 1960.230 1182.960 1972.220 70100 2006.030 1139.650 1966.000 70111 1118.736 1622.800 1640.767 70121 1294.871 1079.156 1618.750 70131 1802.302 1296.800 1963.205 90400 1042.410 664.460 1429.850 -66 1.0 1.0 0.001 0.001 1000.0 .0001 1.0 1.0 1.0 5078 1723.234 1313.129 2692.871 -1 -1 -77 153.48 5079 1791.076 1039.008 2693.495 -1 -1 -76 153.48 5130 1739.959 1311.585 2680.407 1 1 -73 153.48 5131 1831.113 1075.881 2675.282 0 1 -72 153.48 6 70090 70100 1010 1020 1030 1040 7 70090 70100 90400 1010 1020 1030 1040 1010 1729.160 1274.240 1935.650 1020 1583.440 1325.790 1888.610 1030 1556.910 1435.230 1886.820 1040 1631.070 1204.760 1863.880 -66 4.10 950 //GO.FT06F001 DD SYSOUT=* //GO.FT21F001 DD DSN=&&TEMP1,UNIT=DASD,DISP=NEW, // VOL=SER=SYS555,SPACE=(TRK,(20,10),RLSE), // DCB=(RECFM=VBS,LRECL=1920,BLKSIZE=1924) //GO.FT23F001 DD DSN=&&TEMP2,UNIT=DASD,DISP=NEW, // VOL=SER=SYS555,SPACE=(TRK,(20,20),RLSE), // DCB=(RECFM=VBS,LRECL=18700,BLKSIZE=18704) //