C C ************************************************************* C * * C * PROGRAM *GEBAT-CC* * C * * C * FOR CAMERA CALIBRATION & CLOSE-RANGE * C * * C * BY : SABRY EL-HAKIM..1978/80 * C * ---------------------------- * C * * C * * C * * BLOCK-VARIANT..17 PARAMETERS/PHOTO * C * (6 BASIC + 2 P.P. + 1 F + 8 HARMONICS) * C * * COMBINED PHOTO AND DISTANCES ADJUSTMENT * C * * C * LIMITS : * C * -------- * C * - MAX. 4 PHOTOS * C * - MAX. 80 OBJECT POINTS * C * - MAX. 40 GEODETIC POINTS * C * - MAX. 1600 DISTANCES * C * - NO LIMIT ON PHOTO POINTS * C * * C ************************************************************* C C IMPLICIT REAL*8(A-H,O-Z) MAIN INTEGER Q DIMENSION L(321),XP(321),YP(321),LG(321),XG(321),YG(321),ZG(321) *,NF(80),PX2(3,240),LC(080),NV(80),XF(321),YF(321),LQ(080) *,FF(10),PF(10) REAL*8 AP1(2,68),MP(2,2),MPA(2,68),W(2),UP1(68),NP11(17,68),X(68), *M(3,3),BP(2,2),AP2(2,240),NP21(240,68),NP22(3,240),UP2(240), *X1(40,40),X2(40,40),AG(1,68),P(240),NG(68,68),UG(120),XO(120), *NE(68,1),NOG(120,120),VEC(240),R(1,68),B(68),Y(68),PX1(17),XC(80), *YC(80),ZC(80),XQ(80),YQ(80),ZQ(80),AL(321),ALG(321) 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) 47 FORMAT(I10,3F20.9) 409 FORMAT(3F10.5) 489 FORMAT(3F10.5) 201 FORMAT(3F20.3/8F10.1/6F10.4) 410 FORMAT(4I10) 419 FORMAT(3F20.4) 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.****') IRC=5 IWP=6 READ(IRC,410)IL1,IL2,IL3,IJKL READ(IRC,410)IG,IQ READ(IRC,410)NITR,Q IS1=0 IS2=0 IS3=0 IGN=120 CCC=0.00009 2345 FORMAT(5X,16('*')/5X,'*',14X,'*'/5X,'*',' READ-IN DATA ','*'/5X, *'*',14X,'*'/5X,16('*')/) WRITE(IWP,2345) 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 C READING THE PHOTO NO. NF AND POINT NO. L AND PHOTO COORD. XP & YP C READ(IRC,PF)NFIII WRITE(IWP,1234)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)/1000. YP(I)=YP(I)/1000. LOIO=L(I) IF(IG.EQ.0)GO TO 987 ICHEK=L(I)-L(I)/10*10 IF(ICHEK.EQ.IQ)GO TO 985 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 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 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 I=I+1 K=K+1 GO TO 14 13 J=J+1 GO TO 16 15 NP=I-1 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 READ(IRC,1)(PF(I),I=1,10) I=0 55 I=I+1 READ(IRC,PF)LC(I),XC(I),YC(I),ZC(I) IF(LC(I).NE.-66)GO TO 55 IC=I-1 56 FORMAT(I10) DO 762J=1,NP LL=(AL(J)+0.01)/1000000000 A1L=LL AR=AL(J)-A1L*1000000000.0+0.01 IR=AR IR=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) ALG(J)=AL(J) GO TO 762 3102 IF(I.EQ.IC)WRITE(IWP,3103)IR 31 CONTINUE 762 CONTINUE C C READING THE WT. MATRIX OF THE PHOTO POINTS & PX1, PX1U&PX2 C READ(IRC,24)P1,P2,P3 24 FORMAT(3F10.3) PX=P2/(P1*P2-P3**2) PY=P1/(P1*P2-P3**2) PXY=P3/(P1*P2-P3**2) JC=IC NU=17 READ(IRC,201)(PX1(I),I=1,NU) READ(IRC,409)PXC,PYC,PZC READ(IRC,419)PXCO,PYCO,PZCO NUN=NU*NFI MBW=NUN NAN=NUN 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)AOM,PHI,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 I=1,3 DO 3698 J=1,JCS VEC(J)=0.0 3698 PX2(I,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,4034)P1,P2,P3 WRITE(IWP,4035)(PX1(J),J=1,4),(PX1(K),K=12,17) WRITE(IWP,4036)PXC,PYC,PZC WRITE(IWP,4037)PXCO,PYCO,PZCO 4034 FORMAT(/5X,'++WTS. FOR IMAGE COORDS..(PX,PY,PXY):'/10X,3F10.3) 4035 FORMAT(/5X,'++WTS. FOR CAL. & OR. PAR. :'/7X,'XO,YO,F :',3F15.5/ *7X,'HARM. PAR.:',F15.5/7X,'XC,YC,ZC :',3F15.5/7X,'OM,FI,KPA :', *3F15.5) 4036 FORMAT(/5X,'++WTS. FOR OBJECT COORDS. :'/10X,3F15.5) 4037 FORMAT(/5X,'++WTS. FOR CONTROL PTS. :'/10X,3F15.5///) READ(IRC,12)ICC DO 879 I=1,ICC 879 READ(IRC,47)LQ(I),XQ(I),YQ(I),ZQ(I) IDOR=0 1529 DO 524 I=1,JC PX2(1,3*I-2)=PXC PX2(2,3*I-1)=PYC 524 PX2(3,3*I)=PZC 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(2,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(3,3*INV)=PZCO 501 CONTINUE ICHT=0 DO 142 I=1,NU DO 142 J=1,NUN 142 NP11(I,J)=0.0 DO 143 I=1,NAN 143 UP1(I)=0.0 9182 CONTINUE ITR=ITR+1 INB=NP NPS=INB NP=321 INU=NUN INC=JCS JCS=240 NUN=68 NAN=NUN C C CALLING PNORM1 TO FORM NP11, NP21, NP22, UP1 AND UP2 C CALL FNORM1(L,XP,YP,XG,YG,ZG,X,NU,NP,NFI,AP1,NP21,MP,W,JCS, *NP11,INB,MPA,N1,NUN,NAN,UP1,LG,M,BP,PX,PY,PXY,AP2,AL,ALG, *NP22,UP2,MBW,B,VEC,ITR,XF,YF,IDOR,IS2,INC,INU,PX1) NP=INB IF(IDOR.EQ.1)GO TO 1000 IF(IG.EQ.0)GO TO 228 IGM=40 NP=321 CALL SNORM(NOG,UG,NAN,NP,XG,YG,ZG,IG,ITR,X1,X2,INB,IGM,L,XO,IGN, *AG,WG,AL) NP=INB 228 CONTINUE JCS=240 NUN=68 NAN=NUN CALL OORM2(NP22,NOG,PX2,NM,JCS,JC,IG,Y,M,MGW,IGN,3,INC,INU) JCS=240 NUN=68 NAN=NUN CALL OVEC(UP1,NOG,NP21,UG,UP2,NP22,P,JC,IGN,JCS,IG,NNO,VEC,3,INC, *INU) ICR=0 DO 7714 I=1,NFI DO 7714 J=1,NU ICR=ICR+1 7714 UP1(ICR)=UP1(ICR)-B(ICR)*PX1(J) JCS=240 NUN=68 NAN=NUN CALL OORMAL(JC,JCS,NUN,NAN,NG,NP22,NP11,NE,IG,PX1,VEC,NOG,IGN,INU, *NP21,INC,IJKL) JCS=INC NUN=INU NAN=NUN DO 1420 I=1,NAN DO 1420 J=I,NAN 1420 NG(I,J)=NG(J,I) NAN=68 CALL KHLSK(NAN,NG,R,INDEF,ITR,N1,INU,Q) NAN=INU DO 207 I=1,NAN DO 207 J=1,NAN 207 NG(I,J)=0.0 REWINDIL1 DO 205 I=1,NAN READ(IL1)(R(1,II),II=1,NAN) DO 206 J=I,NAN 206 NG(I,J)=R(1,J) 205 CONTINUE NAN=68 CALL KHLSKS(NAN,NG,UP1,B,Y,ITR,N1,INU,Q) NAN=NUN 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) WRITE(IWP,340)I,X(M1),X(M2),X(M3),X(M4),X(M5),X(M6) WRITE(IWP,341)B(M1),B(M2),B(M3),B(M4),B(M5),B(M6) 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(I3,'-ADJ.PAR.=',3F12.3,3F10.4) 341 FORMAT(4X,'RESIDUAL=',3F12.3,3F10.4//) 330 FORMAT(//20X,'XC',10X,'YC',10X,'ZC',7X,'OMEGA',7X,'PHI',5X, *'KAPPA'/) 342 FORMAT(5X,'THE BASIC PAR. + HARMONIC PAR.'/5X,32('-')) DO 40 I=1,JCS 40 P(I)=0.0 DO 402 I=1,NUN DO 402 K=1,JCS 402 P(K)=P(K)+NP21(K,I)*B(I) 1000 CONTINUE DO 709 I=1,JCS 709 UP2(I)=UP2(I)+P(I) JCS=240 NUN=68 NAN=NUN CALL OVEC(UP1,NOG,NP21,UG,UP2,NP22,P,JC,IGN,JCS,IG,NNO,VEC,4,INC, *INU) JCS=INC NUN=INU NAN=NUN WRITE(IWP,20) 20 FORMAT(///5X,'ADJ. COORDS. & CORRN.'/5X,22('*')) 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=(ALG(I)+0.01)/1000000000 A1L=LL AR=ALG(I)-A1L*1000000000.0+0.01 LGOIO=AR LGOIO=LGOIO/1000 LL=(ALG(JJ)+0.01)/1000000000 A1L=LL AR=ALG(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.OR.IJKL.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=321 CALL OHECK(LG,XG,YG,ZG,LQ,XQ,YQ,ZQ,NQ,NP,INB,ICC,ALG) NP=INB IF(IS3.EQ.0)GO TO 5011 IF(ICHT.EQ.0)IDOR=1 IF(ICHT.EQ.0)GO TO 9182 5011 IF(ITR.LT.NITR)GO TO 501 IF(IG.EQ.0)STOP CALL SNORM(NOG,UG,NAN,NP,XG,YG,ZG,IG,0,X1,X2,INB,IGM,L,XO,IGN,AG, *WG,AL) STOP 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 SYINV(A,N1,NN,SE) IMPLICIT REAL*8(A-H,O-Z) REAL*8 A(N1,N1),SE(N1) 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.0D00/EP N=3 200 N=N+1 M=N-1 DO 11 I=1,M SE(I)=0.0D00 DO 10 J=1,M 10 SE(I)=SE(I)+A(I,J)*A(J,N) 11 CONTINUE EP=0.0D00 DO 12 J=1,M 12 EP=EP+A(N,J)*SE(J) EB=A(N,N)-EP EI=1.0D00/EB DO 100 I=1,M DO 100 J=1,M 100 A(I,J)=A(I,J)+SE(I)*SE(J)*EI DO 150 I=1,M A(I,N)=-SE(I)*EI 150 A(N,I)=A(I,N) A(N,N)=EI IF(NN.GT.N)GO TO 200 RETURN END SUBROUTINE OORMAL(JC,JJS,NUM,NAM,NG,NP22,NP11,NE,IG,PX1,VEC,NO, *IG1,INU,NP21,INC,IGL) IMPLICIT REAL*8(A-H,O-Z) REAL*8 NG(NAM,NAM),NP22(3,JJS),NP11(17,NUM),VEC(JJS),PX1(17), *NE(NAM,1),NP21(JJS,NUM),NO(IG1,IG1) JCS=INC NUN=INU NAN=NUN IT=IG+1 NFI=NUN/17 L=1 DO 10 K=1,NAN IF(IG.EQ.0)GO TO 222 DO 20 J=1,IG DO 20 II=1,3 JJ=3*J-3+II VEC(JJ)=0.0 20 CONTINUE DO 15 I=1,IG 15 VEC(JJ)=VEC(JJ)+NP21(3*I-2,K)*NO(JJ,3*I-2)+NP21(3*I-1,K)* *NO(JJ,3*I-1)+NP21(3*I,K)*NO(JJ,3*I-0) 222 CONTINUE DO 25 J=IT,JC DO 25 II=1,3 JJ=3*J-3+II VEC(JJ)=NP21(3*J-2,K)*NP22(II,3*J-2)+NP21(3*J-1,K)*NP22(II,3*J-1) *+NP21(3*J,K)*NP22(II,3*J) 25 CONTINUE DO 40 II=1,NUN NE(II,1)=0.0 DO 41 KK=1,JCS 41 NE(II,1)=NE(II,1)-NP21(KK,II)*VEC(KK) 40 CONTINUE IK=K-(K-1)/17*17 KF=(K-1)/17+1 DO 43 III=1,17 JJ=KF*17-17+III 43 NE(JJ,1)=NP11(III,K)+NE(JJ,1) NE(K,1)=NE(K,1)+PX1(IK) 60 DO 120 I=1,NAN 120 NG(I,K)=NE(I,1) 10 CONTINUE RETURN END SUBROUTINE OHECK(L,X,Y,Z,LC,XC,YC,ZC,IC,NB,INB,IN,AL) IMPLICIT REAL*8 (A-H,O-Z) REAL*8 X(NB),Y(NB),Z(NB),XC(80),YC(80),ZC(80),AL(NB) DIMENSION L(NB),LC(80) IWP=6 NP=INB SX=0. SY=0.0 SZ=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 AD=AA**2 SX=SX+AD BB=Y(K)-B BD=BB**2 SY=SY+BD CC=Z(K)-C 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) WRITE(IWP,21)RMX,RMY,RMZ 21 FORMAT(///5X,'RMS =',3F10.4) 22 FORMAT(I10,3F10.3) RETURN END SUBROUTINE KHLSK(M,A,R,INDEF,ITR,N1,INU,Q) 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 * THE MATRIX R IS STORED ( ROW BY ROW ) IN AN ON LINE DISK * C * * C * MODIFIED BY S. EL-HAKIM AUG9,77 * C ************************************************************************ C IMPLICIT REAL*8(A-H,O-Z) INTEGER Q REAL*8 A(M,M),R(1,M) IWP=6 IL1=21 N=INU REWINDIL1 INDEF=0 DO 30 IP=1,N IF(A(IP,IP))05,05,10 5 WRITE(IWP,6)IP 6 FORMAT(/5X,'ARRAY NOT POSITIVE DEFINITE,COMP#',I4,'IS ADDED TO 1.0 *D-04'/5X,56('$')/) A(IP,IP)=A(IP,IP)+1.0D-4 10 R(01,IP)=DSQRT(A(IP,IP)) IF(IP-N)12,40,40 12 IPPLS1=IP+1 DO 20 K=IPPLS1,N R(01,K)=A(IP,K)/R(01,IP) 20 CONTINUE DO 31 I=IPPLS1,N DO 31 K=I,N 31 A(I,K)=A(I,K)-R(01,I)*R(01,K) WRITE(IL1)(R(1,II),II=1,N) 30 CONTINUE 40 WRITE(IL1)(R(1,II),II=1,N) ENDFILEIL1 RETURN END SUBROUTINE KHLSKS(M,R,B,X,Y,ITR,N1,INU,Q) 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 REAL*8 R(M,M),B(M),X(M),Y(M) N=INU DO 30 K=1,N S=B(K) IF(K-1)25,25,10 10 KLSS1=K-1 DO 20 I=1,KLSS1 S=S+R(I,K)*Y(I) 20 CONTINUE 25 Y(K)=-S/R(K,K) 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 S=S-R(I,K)*X(K) 40 CONTINUE 45 X(I)=S/R(I,I) 50 CONTINUE RETURN END SUBROUTINE OORM2(NP22,NG,PX2,NM,JJS,JC,IG,SE,M,MGW,NAM,IRM,INC *,INU) IMPLICIT REAL*8(A-H,O-Z) REAL*8 NP22(3,JJS),NG(NAM,NAM),PX2(3,JJS),SE(NAM),M(3,3) NM=3*IG JCS=INC NUN=INU NAN=NUN 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(2,III-1) 100 NP22(3,III-0)=NP22(3,III-0)+PX2(3,III) IF(IG.EQ.0)GO TO 1000 DO 200 I=1,IG III=3*I KII=III+2 NG(KII-4,KII-4)=NG(KII-4,KII-4)+NP22(1,III-2) NG(KII-4,KII-3)=NG(KII-4,KII-3)+NP22(1,III-1) NG(KII-4,KII-2)=NG(KII-4,KII-2)+NP22(1,III) NG(KII-3,KII-4)=NG(KII-3,KII-4)+NP22(2,III-2) NG(KII-3,KII-3)=NG(KII-3,KII-3)+NP22(2,III-1) NG(KII-3,KII-2)=NG(KII-3,KII-2)+NP22(2,III) NG(KII-2,KII-4)=NG(KII-2,KII-4)+NP22(3,III-2) NG(KII-2,KII-3)=NG(KII-2,KII-3)+NP22(3,III-1) NG(KII-2,KII-2)=NG(KII-2,KII-2)+NP22(3,III) 200 CONTINUE CALL SYINV(NG,NAM,NM,SE) 1000 CONTINUE IF(IRM.EQ.1.OR.IRM.EQ.4)RETURN K=IG+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 FNORM1(L,XP,YP,XG,YG,ZG,X,NU,NQ,NFI,AP1,NP21,MP,W,JJS *,NP11,INB,MPA,N1,NUM,NAM,UP1,LG,M,BP,PX,PY,PXY,AP2,AL,ALG, *NP22,UP2,MBW,B,VEC,ITR,XF,YF,ICHT,IS2,INC,INU,PX1) C C ************************************************************************ C *SUBROUTINE PNORM1 * C * * C * THIS SUBROUTINE FORMS THE PHOTOGRAMMETRIC NORMAL * C * EQUATIONS ((NP11))-((NP21))-((NP22)) * C * || || || * C * $$$$$$$$$$$$$$$$$$$$$$$$$ * C * AND THE VECTORS UP1 AND UP2 * C * $$$$$$$$$$$ $ C * WHERE THE INPUT JCS=3*NO. OF OBJECT POINTS * C * **CPU = 24 SEC. FOR 30 PHOTOS & 200 POINTS * C * * C * WRITTEN BY S. EL-HAKIM JULY16,77 * C ************************************************************************ C IMPLICIT REAL*8(A-H,O-Z) DIMENSION L(NQ),XP(NQ),YP(NQ),XG(NQ),YG(NQ),ZG(NQ),LG(NQ) REAL*8 AP1(2,NUM),AP2(2,JJS),MP(2,2),MPA(2,NUM),W(2),UP1(NAM), *NP11(17,NUM),XF(NQ),YF(NQ),X(NAM),M(3,3),BP(2,2),NP21(JJS,NAM), *NP22(3,JJS),UP2(JJS),B(NAM),VEC(JJS),PX1(17),AL(NQ),ALG(NQ) IWP=6 IS2=1 WRITE(IWP,50) 50 FORMAT(// 5X,'RESIDUALS AT IMAGE COORDS. :'/5X,28('-')/17X,'PT.', *10X,'VX',13X,'VY'/) JCS=INC NUN=INU NAN=NUN NP=INB NBW=NUN INDX=1 SEGX=0.0 SEGY=0.0 DO 365 I=1,JCS DO 365 J=1,NUN 365 NP21(I,J)=0.0 IF(ITR.GT.1)GO TO 3073 DO 3074 I=1,NP XF(I)=XP(I) 3074 YF(I)=YP(I) 3073 CONTINUE DO1262 J=1,JCS 1262 UP2(J)=0.0 IF(ICHT.EQ.1)GO TO 200 DO 262 J=1,JCS DO 262 I=1,3 262 NP22(I,J)=0.0 IDO=NBW 200 J=1 100 IX=(AL(J)+0.01)/1000000000 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 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 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)) 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=NU-3 NU=17 DVX=T*XR DVY=T*YR IF(ICHT.EQ.1)GO TO 300 AP1(1,KO-5)=-(XR+DVX)*M(3,1)-X(NUOIX-14)*M(1,1) AP1(2,KO-5)=-(YR+DVY)*M(3,1)-X(NUOIX-14)*M(2,1) AP1(1,KO-4)=-(XR+DVX)*M(3,2)-X(NUOIX-14)*M(1,2) AP1(2,KO-4)=-(YR+DVY)*M(3,2)-X(NUOIX-14)*M(2,2) AP1(1,KO-3)=-(XR+DVX)*M(3,3)-X(NUOIX-14)*M(1,3) AP1(2,KO-3)=-(YR+DVY)*M(3,3)-X(NUOIX-14)*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(NUOIX-14)*(-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(NUOIX-14)*(-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(NUOIX-14)*(-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(NUOIX-14)*(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(NUOIX-14)*(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(NUOIX-14)*(-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 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 NUNO3=NUOIX-14 AP2(1,3*KN-2)=(XR+DVX)*M(3,1)+X(NUNO3)*M(1,1) AP2(2,3*KN-2)=(YR+DVY)*M(3,1)+X(NUNO3)*M(2,1) AP2(1,3*KN-1)=(XR+DVX)*M(3,2)+X(NUNO3)*M(1,2) AP2(2,3*KN-1)=(YR+DVY)*M(3,2)+X(NUNO3)*M(2,2) AP2(1,3*KN)=(XR+DVX)*M(3,3)+X(NUNO3)*M(1,3) AP2(2,3*KN)=(YR+DVY)*M(3,3)+X(NUNO3)*M(2,3) NON=NUN NUN=NK+3 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)*(XP(J)*CL+ *R)+X(NUN+11)*XP(J)*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)*XP(J)*SL+ *X(NUN+11)*(R+YP(J)*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 W1=(XR+DVX)*T3+X(NUNO3)*T1 W2=(YR+DVY)*T3+X(NUNO3)*T2 IF(IS2.EQ.0)GO TO 600 WX=W1+AP1(1,KO-5)*B(KO-5)+AP1(1,KO-4)*B(KO-4)+AP1(1,KO-3)*B(KO-3) *+AP1(1,KO-2)*B(KO-2)+AP1(1,KO-1)*B(KO-1)+AP1(1,KO)*B(KO) WY=W2+AP1(2,KO-5)*B(KO-5)+AP1(2,KO-4)*B(KO-4)+AP1(2,KO-3)*B(KO-3) *+AP1(2,KO-2)*B(KO-2)+AP1(2,KO-1)*B(KO-1)+AP1(2,KO)*B(KO) WX=WX-AP2(1,3*KN-2)*VEC(3*KN-2)-AP2(1,3*KN-1)*VEC(3*KN-1) *-AP2(1,3*KN)*VEC(3*KN) WY=WY-AP2(2,3*KN-2)*VEC(3*KN-2)-AP2(2,3*KN-1)*VEC(3*KN-1) *-AP2(2,3*KN)*VEC(3*KN) VX=MP(1,1)*WX+MP(1,2)*WY VY=MP(2,1)*WX+MP(2,2)*WY CX=-PX*(BP(1,1)*VX+BP(2,1)*VY) CY=-PY*(BP(1,2)*VX+BP(2,2)*VY) 2973 FORMAT(I20,2F15.6) WRITE(IWP,2973)ISO,CX,CY SEGX=SEGX+CX**2 SEGY=SEGY+CY**2 600 W(1)=W1*MP(1,1)+W2*MP(1,2) W(2)=W1*MP(2,1)+W2*MP(2,2) IF(ICHT.EQ.1)GO TO 400 DO 110 II=1,2 DO 120 K=1,NU I=IX*NU-NU+K MPA(II,I)=0.0 DO 120 JJ=1,2 120 MPA(II,I)=MPA(II,I)+MP(II,JJ)*AP1(JJ,I) 110 CONTINUE DO 150 II=1,3 I2=KN*3-3+II DO 160 K=1,NU KK=IX*NU-NU+K DO 161 JJ=1,2 161 NP21(I2,KK)=NP21(I2,KK)+AP2(JJ,I2)*MPA(JJ,KK) 160 CONTINUE 150 CONTINUE DO 130 II=1,NU I2=IX*NU-NU+II DO 140 K=1,NU KK=IX*NU-NU+K DO 140 JJ=1,2 140 NP11(II,KK)=NP11(II,KK)+AP1(JJ,I2)*MPA(JJ,KK) 130 CONTINUE DO 210 II=1,NU IK=IX*NU-NU+II DO 210 KK=1,2 210 UP1(IK)=UP1(IK)+AP1(KK,IK)*W(KK) 400 CONTINUE DO 111 II=1,2 DO 121 K=1,3 MPA(II,K)=0.0 KK=KN*3-3+K DO 121 JJ=1,2 121 MPA(II,K)=MPA(II,K)+MP(II,JJ)*AP2(JJ,KK) 111 CONTINUE IF(ICHT.EQ.1)GO TO 500 DO 131 II=1,3 I2=KN*3-3+II DO 141 K=1,3 KK=KN*3-3+K DO 141 JJ=1,2 141 NP22(II,KK)=NP22(II,KK)+AP2(JJ,I2)*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,IK)*W(KK) IF(J.EQ.NP)GO TO 1000 J=J+1 GO TO 100 1000 DF=NP SEGX=DSQRT(SEGX/DF) SEGY=DSQRT(SEGY/DF) WRITE(IWP,1001)SEGX,SEGY 1001 FORMAT(//5X,'ST. ERRORS'/2F20.6//) SEG=SEGX+SEGY DO 1002 I=1,NUN K=I-(I-1)/17*17 1002 SEG=SEG+B(I)**2*PX1(K) DF=NP RETURN END SUBROUTINE OVEC(UP1,NG,NP21,UG,UP2,NP22,Y,JC,NAM,JJS,IG,NNO,VEC,L, *INC,INU) IMPLICIT REAL*8(A-H,O-Z) REAL*8 Y(JJS),UP1(NAM),UP2(JJS),NG(NAM,NAM),UG(NAM),NP22(3,JJS) *,NP21(JJS,NAM),VEC(JJS) 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) DO 10 K=1,IG DO 10 I=1,3 IP=3*K-3+I 10 Y(IP)=Y(IP)+UP2(IP) DO 20 I=1,IG DO 20 K=1,3 JP=3*I-3+K VEC(JP)=0.0 DO 25 J=1,II 25 VEC(JP)=VEC(JP)+NG(JP,J)*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.2.OR.L.EQ.4)RETURN NF=NUN/17 DO 40 II=1,NUN DO 50 K=1,JCS 50 UP1(II)=UP1(II)-NP21(K,II)*VEC(K) 40 CONTINUE RETURN END SUBROUTINE SNORM(NG,UG,NAN,NP,XG,YG,ZG,IG,ITR,S,WT,INB,IGM,L,XO, *IGN,AG,WG,AL) IMPLICIT REAL*8(A-H,O-Z) REAL*8 NG(IGN,IGN),UG(NAN),XG(NP),YG(NP),ZG(NP),S(IGM,IGM), *WT(IGM,IGM),XO(NAN),AG(1,NAN),AL(NP) DIMENSION L(NP) 10 FORMAT(6X,2I5,4X,2F10.3) IRC=5 IWP=6 IF(ITR.EQ.1)WRITE(IWP,11) 11 FORMAT(//5X,'READ IN DIST. OBSERVATIONS'/11X,'FROM',8X,'TO',10X, *'ST.E.',10X,'DIST.'/) 111 FORMAT(//5X,'DISTANCES FROM ADJ. COORDS'/11X,'FROM',8X,'TO',10X, *'DIFF.',10X,'DIST.'/) 12 FORMAT(5X,2I10,5X,F10.3,5X,F10.3) IGS=IG*3 DO 20 II=1,IGS UG(II)=0.0 DO 20 JJ=1,IGS 20 NG(II,JJ)=0.0 IF(ITR.NE.1)GO TO 40 DO 33 K=1,IG DO 33 M=1,IG 33 S(K,M)=0.0 READ(IRC,10)NO DO 77 KM=1,NO READ(IRC,10)I,J,P,SS WRITE(IWP,12)I,J,P,SS IF(I.EQ.-999)GO TO 40 I1=0 I2=0 DO 105 II=1,INB IR=L(II)-L(II)/1000000*1000000 IF(IR.EQ.I)GO TO 104 IF(IR.EQ.J)GO TO 106 GO TO 105 104 LL=(AL(II)+0.01)/1000000000 A1L=LL AR=AL(II)-A1L*1000000000.+.01 IR=AR K=IR-IR/1000*1000 I1=1 KK=3*K XO(KK-2)=XG(II) XO(KK-1)=YG(II) XO(KK)=ZG(II) GO TO 107 106 LL=(AL(II)+0.01)/1000000000 A1L=LL AR=AL(II)-A1L*1000000000.+.01 IR=AR M=IR-IR/1000*1000 I2=1 MM=3*M XO(MM-2)=XG(II) XO(MM-1)=YG(II) XO(MM)=ZG(II) 107 IF(I1.EQ.1.AND.I2.EQ.1)GO TO 108 105 CONTINUE 108 CONTINUE S(K,M)=SS WT(K,M)=1/P*2.0 77 CONTINUE WRITE(IWP,13)IG,NO 13 FORMAT(/5X,'NO. OF GEODETIC PTS. =',I4/5X,'NO. OF GEODETIC OBS. =' *,I4) 40 CONTINUE IF(ITR.EQ.1)GO TO 45 DO 48 II=1,IG III=3*II DO 505 I=1,INB LM=AL(I)/1000. AJ=LM LL=AL(I)-AJ*1000+0.01 IF(LL.EQ.II)GO TO 506 505 CONTINUE 506 CONTINUE XO(III-2)=XG(I) XO(III-1)=YG(I) XO(III-0)=ZG(I) 48 CONTINUE 45 CONTINUE II=0 JJ=0 WRITE(IWP,111) DO 100 KK=1,NO DO 94 I=1,IGS 94 AG(1,I)=0.0 DO 101 I=1,IG DO 101 J=1,IG IF(I.LT.II)GO TO 101 IF(I.LE.II.AND.J.LE.JJ)GO TO 101 IF(I.EQ.J)GO TO 101 IF(S(I,J).EQ.0.0)GO TO 101 I3=3*I J3=3*J SC=DSQRT((XO(I3-2)-XO(J3-2))**2+(XO(I3-1)-XO(J3-1))**2+(XO(I3)- *XO(J3))**2) D=SC-S(I,J) WGOKKO=D*WT(I,J) DO 97 IN=1,INB LM=AL(IN)/1000 LL=AL(IN)-LM*1000+0.01 IF(LL.EQ.I)LOIO=L(IN)-L(IN)/1000000*1000000 IF(LL.EQ.J)LOJO=L(IN)-L(IN)/1000000*1000000 97 CONTINUE WRITE(IWP,12)LOIO,LOJO,D,SC II=I JJ=J IF(ITR.EQ.0)GO TO 100 C1=(XO(J3-2)-XO(I3-2))/SC C2=(XO(J3-1)-XO(I3-1))/SC C3=(XO(J3)-XO(I3))/SC AG(1,I3-2)=-C1*WT(I,J) AG(1,I3-1)=-C2*WT(I,J) AG(1,I3)=-C3*WT(I,J) AG(1,J3-2)=-AG(1,I3-2) AG(1,J3-1)=-AG(1,I3-1) AG(1,J3)=-AG(1,I3) GO TO 701 101 CONTINUE IF(I.EQ.IG.AND.J.EQ.IG)RETURN 701 CONTINUE DO 301 IB=1,3 IA=IB-1 DO 301 JB=1,3 JA=JB-1 I3OJA=I3-JA I3OIA=I3-IA 301 NG(I3OIA,I3OJA)=NG(I3OIA,I3OJA)+AG(1,I3OIA)*AG(1,I3OJA)/WT(I,J) DO 302 IB=1,3 IA=IB-1 DO 302 JB=1,3 JA=JB-1 J3OJA=J3-JA I3OIA=I3-IA 302 NG(I3OIA,J3OJA)=NG(I3OIA,J3OJA)+AG(1,I3OIA)*AG(1,J3OJA)/WT(I,J) DO 303 IB=1,3 IA=IB-1 DO 303 JB=1,3 JA=JB-1 J3OIA=J3-IA I3OJA=I3-JA 303 NG(J3OIA,I3OJA)=NG(I3OJA,J3OIA) DO 304 IB=1,3 DO 304 JB=1,3 IA=IB-1 JA=JB-1 J3OJA=J3-JA J3OIA=J3-IA 304 NG(J3OIA,J3OJA)=NG(J3OIA,J3OJA)+AG(1,J3OIA)*AG(1,J3OJA)/WT(I,J) DO 305 IB=1,3 IA=IB-1 305 UG(I3-IA)=UG(I3-IA)+AG(1,I3-IA)*WGOKKO/WT(I,J) DO 306 IB=1,3 IA=IB-1 306 UG(J3-IA)=UG(J3-IA)+AG(1,J3-IA)*WGOKKO/WT(I,J) 100 CONTINUE RETURN END