//PTBV JOB /*SERVICE -4 /*JOBPARM S=20,L=15,R=4096,PRINT=ALL // EXEC FORTVCLG,RC=4096K,RL=1024K,RG=1024K,LIB=LIBD, // PARM.FORT='NOMAP,NOXREF' //FORT.SYSIN DD * C C *************************************************************** C * PROGRAM PTBV * C * * C * 'PHOTOGRAMMETRIC TRIANGULATION BY BUNDLES-PHOTO VARIANT' * C * * C * WRITTEN BY: COSTAS ARMENAKIS ON JULY 1986 * C *************************************************************** C IMPLICIT REAL*8 (A-H,O-Z) C C IR = INPUT UNIT C IW = OUTPUT UNIT C IEX1 = EXTERNAL UNIT 1 C IEX2 = EXTERNAL UNIT 2 C IBYT1 = NUMBER OF BYTES PER RECORD C IBYT2 = NUMBER OF BYTES PER RECORD C MAXPH = MAXIMUM NUMBER OF PHOTOGRAPHS C MAXPPP = MAXIMUM NUMBER OF POINTS PER PHOTOGRAPH C NAPL = TOTAL NUMBER OF PHOTO-POINTS C MAXOP = MAXIMUM NUMBER OF OBJECT POINTS C MAXPA = MAXIMUN NUMBER OF PHOTOGRAPHS ON WHICH A PHOTO-POINT C APPEARS PLUS 3 C NPHOT = TOTAL NUMBER OF PHOTOGRAPHS C JFP = TOTAL NUMBER OF OBJECT POINTS C NGF = TOTAL NUMBER OF CONTROL POINTS (H&V) C NC = TOTAL NUMBER OF CHECK POINTS C XG,YG,ZG= CONTROL POINTS C XC,YC,ZC= CHECK POINTS C ITRF = TOTAL NUMBER OF ITERATIONS C ITRI = NUMBER OF ITERATION TO INTRODUCE THE IO PARAMETERS C COMMON /DEV/ IR,IW,IEX1,IEX2,IBYT1,IBYT2 COMMON /DIM/ MAXPH,MAXPPP,NALP,MAXOP,MAXPA COMMON /ARR/ XL(50),YL(50),DID(100),XG(100),YG(100),ZG(100), @ WX(100),WY(100),WZ(100),XK(100),YK(100),ZK(100), @ IPT(50),IG(100),IK(100), @ INDEX(100),IAPN(450),IIPH(450),LTAB(1000,12) COMMON /VAR/ PX,PY,PAUX,NPHOT,JFP,NGF,ITRF,ITRI,ITER,NC,LDEF COMMON /ERR/ IERR1,IERR2,IERR3,IERR4 C IR=5 IW=6 IEX1=20 IEX2=21 IBYT1=2200 IBYT2=80 MAXPH=9 MAXPPP=50 NALP=450 MAXOP=1000 MAXPA=12 C OPEN (UNIT=IEX1,ERR=9000,STATUS='NEW',FILE='PHOTDAT', @ ACCESS='DIRECT',FORM='UNFORMATTED',RECL=IBYT1) C OPEN (UNIT=IEX2,ERR=9002,STATUS='NEW',FILE='APRXOBJ', @ ACCESS='DIRECT',FORM='UNFORMATTED',RECL=IBYT2) C CALL PHOTO IF (IERR1.EQ.1) GO TO 9900 C CALL OBJECT IF (IERR2.EQ.2) GO TO 9900 IF (IERR3.EQ.3) GO TO 9900 IF (IERR4.EQ.4) GO TO 9900 C CALL ADJUST C GO TO 9900 9000 WRITE(IW,9001) 9001 FORMAT('0',10X,'ERROR DURING OPENING OF FILE ',/) GO TO 9900 9002 WRITE(IW,9003) 9003 FORMAT('0',10X,'ERROR DURING OPENING OF FILE ',/) 9900 CLOSE (UNIT=IEX1,ERR=9999,STATUS='DELETE') C - - CLOSE (UNIT=IEX2,ERR=9999,STATUS='DELETE') CLOSE (UNIT=IEX2,ERR=9999,STATUS='KEEP') 9999 STOP END C C *************************************************************** SUBROUTINE PHOTO C * * C * PROCESS PHOTO-INFORMATION - COMPUTES APPROX. OBJECT COORDS * C * * C * WRITTEN BY: COSTAS ARMENAKIS ON JULY 1986 * C *************************************************************** C IMPLICIT REAL*8 (A-H,O-Z) C COMMON /DEV/ IR,IW,IEX1,IEX2,IBYT1,IBYT2 COMMON /DIM/ MAXPH,MAXPPP,NALP,MAXOP,MAXPA COMMON /ARR/ XL(50),YL(50),DID(100),XG(100),YG(100),ZG(100), @ WX(100),WY(100),WZ(100),XK(100),YK(100),ZK(100), @ IPT(50),IG(100),IK(100), @ INDEX(100),IAPN(450),IIPH(450),LTAB(1000,12) COMMON /CAM/ OE(6),OI(10),OIUD(10),WOI(10),R(3,3) COMMON /VAR/ PX,PY,PAUX,NPHOT,JFP,NGF,ITRF,ITRI,ITER,NC,LDEF COMMON /ERR/ IERR1,IERR2,IERR3,IERR4 C C PHOTO-SEQUENCE C JPH=1 C C READ IN PHOTO NUMBER C 50 READ(IR,*) IPH IF (IPH.LT.9999990) THEN WRITE(IW,*) IPH C C READ IN ELEMENTS OF INTERIOR ORIENTATION, IN AND THEN XO,YO,C==> C READ(IR,*) (OI(K),K=1,10) DO 60 IL=1,3 OI(IL)=OI(IL)/1000.D0 60 CONTINUE WRITE(IW,*) (OI(K),K=1,10) C COPY TO WHICH CONTAINS THE UPDATED DO 55 IL=1,10 OIUD(IL)=OI(IL) 55 CONTINUE C C READ IN WEIGHTS RELATED TO THE ELEMENTS OF IO C READ(IR,*) (WOI(K),K=1,10) WRITE(IW,*) (WOI(K),K=1,10) C C READ IN ELEMENTS OF EXTERIOR ORIENTATION C POSITION IN / ATTITUDE IN C READ(IR,*) (OE(K),K=1,6) C C TRANSFER ATTITUDE INTO C RAD=DARCOS(-1.D0)/180.D0 DO 70 K=4,6 OE(K)=OE(K)*RAD 70 CONTINUE WRITE(IW,*) (OE(K),K=1,6) C C READ IN PNT ID # AND PHOTO-OBSERVATIONS, = COUNTER C IPN=1 100 READ(IR,*) IPT(IPN),XL(IPN),YL(IPN) IF (IPT(IPN).GT.0) THEN XL(IPN)=XL(IPN)/1000.D0 YL(IPN)=YL(IPN)/1000.D0 WRITE(IW,*) IPT(IPN),XL(IPN),YL(IPN) IPN=IPN+1 IF(IPN.EQ.51) THEN WRITE(IW,8888) IPH IERR1=1 RETURN ENDIF C C READ NEXT POINT C GO TO 100 ELSE C C END OF DATA FOR EACH PHOTO C SORT IMAGE POINTS FOR EACH PHOTO C IPN=IPN-1 DO 120 J=1,IPN DID(J)=DFLOAT(IPT(J)) 120 CONTINUE C CALL SORTD (DID,IPN,1,IPN,1,INDEX,2,XL,YL) C DO 130 J=1,IPN IPT(J)=IDINT(DID(J)) 130 CONTINUE C C CHECK FOR DUBLICATE PNTS / PERFORM LINEAR SEARCH C SET ZERO THE PNT ID # IF PNT APPEARS MORE THAN ONCE AND SEND C WARNING MESSAGE THAT ONLY THE FIRST PNT IS USED C DO 160 L=1,IPN IDUB=IPT(L) IC=0 DO 140 J=1,IPN IF (IPT(J).EQ.IDUB) THEN IC=IC+1 IF (IC.GT.1) THEN MQ=IPT(J) IPT(J)=0 WRITE(IW,5000) MQ,IPH ENDIF ENDIF IF (IPT(J).GT.IDUB) GO TO 160 140 CONTINUE 160 CONTINUE C C ELIMINATE DUBLE POINTS/PHOTO,, = FINAL NUMBER OF PNTS/PHOTO C IPNF=0 DO 180 J=1,IPN IF (IPT(J).NE.0) THEN IPNF=IPNF+1 IPT(IPNF)=IPT(J) XL(IPNF)=XL(J) YL(IPNF)=YL(J) ENDIF 180 CONTINUE C C STORE ALL DATA RELATED TO EACH PHOTO ON ONE RECORD C WRITE (UNIT=IEX1,REC=JPH,ERR=9900) IPH, @ JPH, @ (OIUD(K),K=1,10), @ (OI(K),K=1,10), @ (WOI(K),K=1,10), @ (OE(K),K=1,6), @ IPNF, @ (IPT(K),XL(K),YL(K),K=1,IPNF) WRITE (IW,*) IPH, @ JPH, @ (OIUD(K),K=1,10), @ (OI(K),K=1,10), @ (WOI(K),K=1,10), @ (OE(K),K=1,6), @ IPNF, @ (IPT(K),XL(K),YL(K),K=1,IPNF) C C GO FOR THE NEXT PHOTO C JPH=JPH+1 GO TO 50 ENDIF C C CHECK FOR MINIMUN # OF PHOTOS C ELSE MINP=JPH-1 IF (MINP.LT.2) THEN WRITE(IW,5010) MINP IERR1=1 RETURN ENDIF ENDIF C C INPUT WEIGHTS FOR PHOTO-COORDINATES READ(IR,*) PX,PY C C DETERMINE THE NUMBER OF ALL PHOTO PNTS AND THE PHOTO THEY APPEAR ON C CREATE VECTOR CONTAINING ALL PNT ID NUMBERS C = VECTOR CONTAINING ALL PHOTO PNT ID # ==> EACH ELEMENT C APPEARS MORE THAN TWO TIMES C = VECTOR CONTAINING THE SEQUENTIAL PHOTO-NUMBER WITH WHICH C THE POINT IS RELATED C NPHOT = NUMBER OF PHOTOGRAPHS C IDIM=0 KK=1 MM=0 NPHOT=JPH-1 DO 220 J=1,NPHOT READ(UNIT=IEX1,REC=J,ERR=9900) IPH, @ JPH, @ (OIUD(K),K=1,10), @ (OI(K),K=1,10), @ (WOI(K),K=1,10), @ (OE(K),K=1,6), @ IPNF, @ (IPT(K),XL(K),YL(K),K=1,IPNF) MM=MM+IPNF I=0 DO 200 L=KK,MM I=I+1 IDIM=IDIM+1 IAPN(L)=IPT(I) IIPH(L)=JPH 200 CONTINUE KK=MM+1 220 CONTINUE C C DETERMINE WHICH POINT(S) APPEAR ON WHICH PHOTOGRAPH IN TABULAR C FORM. MAX NUMBER OF PHOTOS ON WHICH A POINT APPEARS IS 9 ==> C DIM=(# OF PNTS,12). FOR : 1ST COL=PNT ID #, C 2ND COL=SEQUENTIAL # OF PHOTOS, 11TH COL=TOTAL # OF PNT APPEARENCE C 12TH COL=CODE # FOR THE TYPE OF CONTROL C = NUMBER OF TIMES BY WHICH EACH PNTS APPEARS C = TOTAL NUMBER OF OBJECT POINTS C JFP=0 DO 280 J=1,IDIM IF (IAPN(J).NE.0) THEN IELEM=IAPN(J) ICO=0 DO 260 K=1,IDIM IF (IAPN(K).EQ.IELEM) THEN LT1=IAPN(K) IAPN(K)=0 ICO=ICO+1 IF (ICO.EQ.1) THEN JFP=JFP+1 LTAB(JFP,1)=LT1 LTAB(JFP,2)=IIPH(K) ELSE LTAB(JFP,ICO+1)=IIPH(K) ENDIF LTAB(JFP,11)=ICO ENDIF 260 CONTINUE ENDIF 280 CONTINUE DO 989 I=1,JFP WRITE(IW,*) (LTAB(I,J),J=1,12) 989 CONTINUE C C DETERMINE APPROXIMATE OBJECT COORDINATES FOR EACH PHOTO-POINT C USE THE FIRST TWO PHOTOS ON WHICH THE POINT APPEARS C FOR EACH POINT IN DO C DO 300 I=1,JFP READ (UNIT=IEX1,REC=LTAB(I,2),ERR=9900) IPH, @ JPH, @ (OIUD(K),K=1,10), @ (OI(K),K=1,10), @ (WOI(K),K=1,10), @ (OE(K),K=1,6), @ IPNF, @ (IPT(K),XL(K),YL(K),K=1,IPNF) CALL FIND (IPNF,IPT,LTAB(I,1),JJ) CALL ROTMAT (OE(4),OE(5),OE(6),SOM,SPH,SKP,COM,CPH,CKP,R) C XCAM1=OE(1) YCAM1=OE(2) ZCAM1=OE(3) U1L=XL(JJ)-OI(1) V1L=YL(JJ)-OI(2) U1=R(1,1)*U1L+R(2,1)*V1L-R(3,1)*OI(3) V1=R(1,2)*U1L+R(2,2)*V1L-R(3,2)*OI(3) W1=R(1,3)*U1L+R(2,3)*V1L-R(3,3)*OI(3) READ (UNIT=IEX1,REC=LTAB(I,3),ERR=9900) IPH, @ JPH, @ (OIUD(K),K=1,10), @ (OI(K),K=1,10), @ (WOI(K),K=1,10), @ (OE(K),K=1,6), @ IPNF, @ (IPT(K),XL(K),YL(K),K=1,IPNF) CALL FIND (IPNF,IPT,LTAB(I,1),JJ) CALL ROTMAT (OE(4),OE(5),OE(6),SOM,SPH,SKP,COM,CPH,CKP,R) C U2L=XL(JJ)-OI(1) V2L=YL(JJ)-OI(2) U2=R(1,1)*U2L+R(2,1)*V2L-R(3,1)*OI(3) V2=R(1,2)*U2L+R(2,2)*V2L-R(3,2)*OI(3) W2=R(1,3)*U2L+R(2,3)*V2L-R(3,3)*OI(3) BX=OE(1)-XCAM1 BY=OE(2)-YCAM1 BZ=OE(3)-ZCAM1 C COMPONENTS OF THE RESIDUAL PARALLAX VECTOR,EX,EY,EZ EX=V1*W2-W1*V2 EY=W1*U2-U1*W2 EZ=U1*V2-V1*U2 C SOLVE THE LINEAR SYSTEM ==> SCALAR MULTIPLIERS DTR=U1*(EZ*V2-EY*W2) @ -EX*(W1*V2-V1*W2) @ -U2*(V1*EZ-W1*EY) DTR1=BX*(EZ*V2-EY*W2) @ -EX*(BZ*V2-BY*W2) @ -U2*(BY*EZ-BZ*EY) DTR2=U1*(BZ*V2-BY*W2) @ -BX*(W1*V2-V1*W2) @ -U2*(V1*BZ-W1*BY) C SCALAR MULTIPLIERS SL1=DTR1/DTR D=DTR2/DTR C COMPUTE APPROXIMATE OBJECT COORDINATES XJ=XCAM1+SL1*U1+0.5*D*EX YJ=YCAM1+SL1*V1+0.5*D*EY ZJ=ZCAM1+SL1*W1+0.5*D*EZ C COMPUTE OBJECT RESIDUAL PARALLAXES RP=D*DSQRT(EX**2+EY**2+EZ**2) WRITE(IW,*) LTAB(I,1),XJ,YJ,ZJ,RP WRITE(UNIT=IEX2,REC=I,ERR=9900) LTAB(I,1),XJ,YJ,ZJ 300 CONTINUE 5000 FORMAT('0',10X,'PHOTO-POINT # ',I6,2X,'ON PHOTOGRAPH ',I6, @ 2X,'APPEARS MORE THAN ONCE. ONLY THE 1ST POINT USED',/) 5010 FORMAT('0',10X,'MORE THAN ',I4,2X,'PHOTOGRAPHS ARE REQUIRED',/) 8888 FORMAT(' ',/,15X,'MORE THAN 50 POINTS IN PHOTO # ',I10) 9900 RETURN END C C *************************************************************** SUBROUTINE ROTMAT (OM,PH,PK,SOM,SPH,SKP,COM,CPH,CKP,R) C * * C * FORMS THE ROTATION MATRIX (R ROTATES THE REFERENCE SYSTEM) * C * * C * WRITTEN BY: COSTAS ARMENAKIS ON JULY 1986 * C *************************************************************** C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION R(3,3) C SOM=DSIN(OM) SPH=DSIN(PH) SKP=DSIN(PK) COM=DCOS(OM) CPH=DCOS(PH) CKP=DCOS(PK) C R(1,1)=CPH*CKP R(1,2)=COM*SKP+SOM*SPH*CKP R(1,3)=SOM*SKP-COM*SPH*CKP R(2,1)=-CPH*SKP R(2,2)=COM*CKP-SOM*SPH*SKP R(2,3)=SOM*CKP+COM*SPH*SKP R(3,1)=SPH R(3,2)=-SOM*CPH R(3,3)=COM*CPH C RETURN END C C ******************************************************************* SUBROUTINE FIND (M,MVEC,NUM,LOC) C * * C * FINDS LOCATION OF IN VECTOR * C * WRITTEN BY: COSTAS ARMENAKIS ON JULY 1986 * C ******************************************************************* C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION MVEC(100) LOC=0 DO 10 L=1,M IF (MVEC(L).EQ.NUM) THEN LOC=L RETURN ENDIF 10 CONTINUE RETURN END C C**************************************************************** SUBROUTINE OBJECT C * * C * PROCESS OBJECT INFORMATION - PERFORMS NECESSARY CHECKS * C * COMPUTES DEGREES OF FREEDOM * C * WRITTEN BY: COSTAS ARMENAKIS ON JULY 1986 * C**************************************************************** C IMPLICIT REAL*8 (A-H,O-Z) C COMMON /DEV/ IR,IW,IEX1,IEX2,IBYT1,IBYT2 COMMON /DIM/ MAXPH,MAXPPP,NALP,MAXOP,MAXPA COMMON /ARR/ XL(50),YL(50),DID(100),XG(100),YG(100),ZG(100), @ WX(100),WY(100),WZ(100),XK(100),YK(100),ZK(100), @ IPT(50),IG(100),IK(100), @ INDEX(100),IAPN(450),IIPH(450),LTAB(1000,12) COMMON /VAR/ PX,PY,PAUX,NPHOT,JFP,NGF,ITRF,ITRI,ITER,NC,LDEF COMMON /HV/ JF,JH,JV COMMON /ERR/ IERR1,IERR2,IERR3,IERR4 C C WEIGHT FACTOR FOR THE APPROXIMATE OBJECT COORDS READ(IR,*) PAUX,ITRF,ITRI WRITE(IW,*) ' PAUX = ',PAUX C C INPUT OBJECT CONTROL POINTS WITH THEIR ID # AND WEIGHTS. IF C A COORDINATE AND THE CORRESPONDING WEIGHT ARE ZERO THE POINT C IS TREATED AS HORIZONTAL OR VERTICAL C NG=1 20 READ(IR,*) IG(NG),XG(NG),YG(NG),ZG(NG),WX(NG),WY(NG),WZ(NG) IF (IG(NG).GT.0) THEN NG=NG+1 GO TO 20 ELSE NG=NG-1 C C SORT THE OBJECT CONTROL POINTS C DO 30 I=1,NG DID(I)=DFLOAT(IG(I)) 30 CONTINUE C CALL SORTD (DID,NG,1,NG,1,INDEX,6,XG,YG,ZG,WX,WY,WZ) C DO 40 I=1,NG IG(I)=IDINT(DID(I)) WRITE(IW,*) IG(I),XG(I),YG(I),ZG(I),WX(I),WY(I),WZ(I) 40 CONTINUE ENDIF C C CHECK FOR DUBLICATE PNTS / PERFORM LINEAR SEARCH C SET ZERO THE PNT ID # IF PNT APPEARS MORE THAN ONCE AND SEND C WARNING MESSAGE THAT ONLY THE FIRST PNT IS USED C DO 80 L=1,NG IDUB=IG(L) IC=0 DO 60 J=1,NG IF (IG(J).EQ.IDUB) THEN IC=IC+1 IF (IC.GT.1) THEN LR=IG(J) IG(J)=0 WRITE(IW,5000) LR ENDIF ENDIF IF (IG(J).GT.IDUB) GO TO 80 60 CONTINUE 80 CONTINUE C ELIMINATE DUBLE CONTRL POINTS, = FINAL NUMBER OF CONTRL PNTS NGF=0 DO 100 J=1,NG IF (IG(J).NE.0) THEN NGF=NGF+1 IG(NGF)=IG(J) XG(NGF)=XG(J) YG(NGF)=YG(J) ZG(NGF)=ZG(J) WX(NGF)=WX(J) WY(NGF)=WY(J) WZ(NGF)=WZ(J) ENDIF 100 CONTINUE C C CHECK THAT ALL CONTROL POINTS APPEAR AS PHOTO-POINTS C ENTER CODE # BASED ON THE TYPE OF CONTROL POINT (1=F,2=H,3=V) C LCOM=0 JF=0 JH=0 JV=0 DO 120 I=1,NGF DO 110 J=1,JFP IF (IG(I).EQ.LTAB(J,1)) THEN LCOM=LCOM+1 IF (XG(I).NE.0.0) THEN IF (ZG(I).NE.0.0) THEN JF=JF+1 LTAB(J,12)=1 ELSE JH=JH+1 LTAB(J,12)=2 ENDIF ELSE JV=JV+1 LTAB(J,12)=3 ENDIF GO TO 120 ENDIF 110 CONTINUE 120 CONTINUE C NOBP = NUMBER OF PHOTOGRAMMETRIC OBS NOBP=0 DO 125 I=1,JFP NOBP=NOBP+2*LTAB(I,11) WRITE(IW,*) (LTAB(I,J),J=1,12) 125 CONTINUE C CHECK IF ALL CONTROL POINTS APPEAR AS PHOTO-POINTS C IF (LCOM.NE.NGF) THEN WRITE(IW,5010) IERR2=2 RETURN ENDIF C C CHECK THE MINIMUN NUMBER OF REFERENCE POINTS REQUIRED TO C DETERMINE THE DATUM C IF (JF.LT.3) THEN IF (JH.LT.2.OR.JV.LT.3) THEN WRITE(IW,5020) IERR3=3 RETURN ENDIF ENDIF C C CHECK IF THERE ANY SINGLE PHOTO-POINTS IN . C IF YES CHECK IF THE CORRESPONDING CONTROL POINTS EXIST. C DO 140 I=1,JFP IF (LTAB(I,11).EQ.1) THEN ISNGL=LTAB(I,1) DO 130 J=1,NGF IF (ISNGL.EQ.IG(J)) THEN LOK=1 ENDIF 130 CONTINUE IF (LOK.NE.1) THEN WRITE(IW,5030) LTAB(I,1) IERR4=4 RETURN ENDIF ENDIF 140 CONTINUE C NOBC = NUMBER OF CONTROL PNT OBS NOBC=0 C REPLACE APPROXIMATE COORDS OF CONTROL PNTS WITH THEIR REAL VALUES C DO 144 IL=1,NGF DO 142 JR=1,JFP IF (IG(IL).EQ.LTAB(JR,1)) THEN READ (UNIT=IEX2,REC=JR,ERR=4000) IDO,XJ,YJ,ZJ IF (LTAB(JR,12).EQ.1) THEN NOBC=NOBC+3 XJ=XG(IL) YJ=YG(IL) ZJ=ZG(IL) ELSEIF (LTAB(JR,12).EQ.2) THEN NOBC=NOBC+2 XJ=XG(IL) YJ=YG(IL) ELSEIF (LTAB(JR,12).EQ.3) THEN NOBC=NOBC+1 ZJ=ZG(IL) ENDIF WRITE(UNIT=IEX2,REC=JR,ERR=4000) IDO,XJ,YJ,ZJ GO TO 144 ENDIF 142 CONTINUE 144 CONTINUE C CHECK IF ALL OBJECT PNTS ARE USED AS OBS IF (PAUX.NE.0.D0) THEN NOBC=3*JFP ENDIF C CHECK IF THE IO PARAMETERS ARE USED AS OBS NOBI=0 IF (ITRI.LE.ITRF) THEN NOBI=10*NPHOT ENDIF C NOB = NUMBER OF OBS NOB=NOBP+NOBC+NOBI C LUNP = NUMBER OF UNKNOWN PARAMETERS LUNP=6*NPHOT+NOBI+3*JFP C LDEF = DEGREES OF FREEDOM LDEF=NOB-LUNP WRITE(IW,*) ' NOB=',NOB,' LUNP=',LUNP,' LDEF=',LDEF C INPUT CHECK POINTS. = # OF CHECK POINTS NC=1 150 READ(IR,*) IK(NC),XK(NC),YK(NC),ZK(NC) IF (IK(NC).GT.0) THEN NC=NC+1 GO TO 150 ELSE NC=NC-1 C C SORT THE CHECK POINTS C DO 160 I=1,NC DID(I)=DFLOAT(IK(I)) 160 CONTINUE C CALL SORTD (DID,NC,1,NC,1,INDEX,3,XK,YK,ZK) C DO 180 I=1,NC IK(I)=IDINT(DID(I)) WRITE(IW,*) IK(I),XK(I),YK(I),ZK(I) 180 CONTINUE ENDIF RETURN 4000 WRITE(IW,4001) 4001 FORMAT('0',10X,'PROBLEMS WITH FILE ',/) 5000 FORMAT('0',10X,'CONTROL POINT # ',I6,3X,'APPEARS MORE THAN ONCE') 5010 FORMAT('0',10X,'ALL CONTROL POINTS DO NOT APPEAR AS PHOTO-POINTS') 5020 FORMAT('0',10X,'NOT ENOUGH CONTROL POINTS') 5030 FORMAT('0',10X,'PHOTO-POINT # ',I6,3X,'IS A SINGLE POINT') RETURN END C C *************************************************************** SUBROUTINE ADJUST C * * C * PERFORMS THE BUNDLE-BLOCK ADJUSTMENT * C * WRITTEN BY: COSTAS ARMENAKIS ON JULY 1986 * C *************************************************************** C IMPLICIT REAL*8 (A-H,O-Z) C INTEGER*4 IPTI(50) C REAL*8 E(2,6),A(2,3),P(2,10),AA(2,3),C(3,3),AAT(3,2), @ C1(3,3),W(3),ET(6,2),E1(6,6),E2(6,3),E3(6,3),E4(6,6), @ E2T(3,6),E0(6,6),PT(10,2),P1(10,6),P2(10,3),P3(10,3), @ AT(3,2),P4(3,6),P5(10,6),P0(10,6),Q0(10,10),Q(10,10), @ Q1(3,10),Q2(10,10),AWOI(10),N(10440) C @ NEO(465),NIO(1275) C REAL*8 WP(2),WPA(2),WC(3),AW1(3),AW(3),EW(6),AW0(6),AW2(6), @ WI(10),QW(10),QW0(10),QW1(10),U(144), @ AOIUD(10),AOI(10),AOE(6),XLI(50),YLI(50),R1(3,3) C @ UEO(54),UIO(90) C COMMON /DEV/ IR,IW,IEX1,IEX2,IBYT1,IBYT2 COMMON /DIM/ MAXPH,MAXPPP,NALP,MAXOP,MAXPA COMMON /ARR/ XL(50),YL(50),DID(100),XG(100),YG(100),ZG(100), @ WX(100),WY(100),WZ(100),XK(100),YK(100),ZK(100), @ IPT(50),IG(100),IK(100), @ INDEX(100),IAPN(450),IIPH(450),LTAB(1000,12) COMMON /CAM/ OE(6),OI(10),OIUD(10),WOI(10),R(3,3) COMMON /VAR/ PX,PY,PAUX,NPHOT,JFP,NGF,ITRF,ITRI,ITER,NC,LDEF COMMON /ERR/ IERR1,IERR2,IERR3,IERR4 C C DIMENSIONS OF MATRICES , , , , AND C JDMN1=6*MAXPH C JDMN2=10*MAXPH C IDMN1=(JDMN1*(JDMN1+1))/2 C IDMN2=(JDMN2*(JDMN2+1))/2 NUDIM=16*MAXPH NDIM=(NUDIM*(NUDIM+1))/2 C INITIALIZE TO ZERO , C DO 6005 JR=1,465 C NEO(JR)=0.D0 C NIO(JR)=0.D0 C6005 CONTINUE C DO 6006 IL=466,1275 C NIO(IL)=0.D0 C6006 CONTINUE C DO 6007 JR=1,30 C UEO(JR)=0.D0 C UIO(JR)=0.D0 C6007 CONTINUE C DO 6008 IL=31,50 C UIO(IL)=0.D0 C6008 CONTINUE C DO 6000 IL=1,NUDIM U(IL)=0.D0 N(IL)=0.D0 6000 CONTINUE DO 6001 IL=NUDIM+1,NDIM N(IL)=0.D0 6001 CONTINUE C ITER=1 6010 CONTINUE WRITE(IW,*) ' ******** ITERATION # ',ITER,' ********' C WRITE(IW,*) 'NPHOT=',NPHOT DO 3000 I=1,NPHOT C DO 2000 J=1,NPHOT C INITIALIZE TO ZERO ,,,, DO 2 JR=1,6 AW0(JR)=0.D0 DO 1 IL=1,6 E0(IL,JR)=0.D0 1 CONTINUE 2 CONTINUE DO 4 JR=1,6 DO 3 IL=1,10 P0(IL,JR)=0.D0 3 CONTINUE 4 CONTINUE DO 6 JR=1,10 QW0(JR)=0.D0 DO 5 IL=1,10 Q0(IL,JR)=0.D0 5 CONTINUE 6 CONTINUE C IF (J.EQ.I) THEN C READ (UNIT=IEX1,REC=I,ERR=9900) IPH, @ JPH, @ (OIUD(IL),IL=1,10), @ (OI(IL),IL=1,10), @ (WOI(IL),IL=1,10), @ (OE(IL),IL=1,6), @ IPNF, @ (IPT(IL),XL(IL),YL(IL),IL=1,IPNF) CALL ROTMAT (OE(4),OE(5),OE(6),SOM,SPH,SKP,COM,CPH,CKP,R) C CFORM DIAGONAL ELEMENTS OF NORMAL EQUATIONS C COMPUTE PHOTO-PNT CONTRIBUTIONS TO /PHOTO KKK=IPNF DO 1200 K=1,KKK DO 10 L=1,JFP READ(UNIT=IEX2,REC=L,ERR=3500) IDO,XJ,YJ,ZJ IF (IPT(K).EQ.IDO) THEN DX=XJ-OE(1) DY=YJ-OE(2) DZ=ZJ-OE(3) GO TO 12 ENDIF 10 CONTINUE 12 CONTINUE U1=R(1,1)*DX+R(1,2)*DY+R(1,3)*DZ U2=R(2,1)*DX+R(2,2)*DY+R(2,3)*DZ U3=R(3,1)*DX+R(3,2)*DY+R(3,3)*DZ T=U3**2 C FORM MATRIX (PART.DERIV. WRT EO ELEMENTS) CALL EMAT (OIUD(3),U1,U2,U3,T,DX,DY,DZ,R,SOM,SPH,SKP, @ COM,CPH,CKP,E) C C FORM MATRIX (PART.DERIV. WRT OBJECT COORDS) CALL AMAT (OIUD(3),U1,U2,U3,T,DX,DY,DZ,R,A) C C FORM

MATRIX (PART.DERIV. WRT IO ELEMENTS) CALL PMAT (OIUD,U1,U2,U3,XL(K),YL(K),P) C FORM MATRIX (MISCLOSURE VECTOR OF PHOTOGRAMMETRIC OBSERVATIONS) CALL WPMAT (OIUD,U1,U2,U3,XL(K),YL(K),WP) C FORM MATRIX (MISCLOSURE VECTOR OF INTERIOR ORIENT OBSERVATIONS) CALL WIMAT (OIUD,OI,WI) C C FORM MATRIX FOR THE KTH PNT. USE MATRIX INSTEAD OF C LOCATE PNT K IN C DO 20 M=1,JFP IF (IPT(K).EQ.LTAB(M,1)) GO TO 30 20 CONTINUE 30 CONTINUE C INITIALIAZE AND FOR THE KTH PNT DO 60 JR=1,3 AW(JR)=0.D0 DO 40 IL=1,3 C(IL,JR)=0.D0 40 CONTINUE 60 CONTINUE C FIND PHOTOS ON WHICH PNT K APPEARS DO 120 IP=1,LTAB(M,11) READ (UNIT=IEX1,REC=LTAB(M,IP+1),ERR=9900) IPHI, @ JPHI, @ (AOIUD(IL),IL=1,10), @ (AOI(IL),IL=1,10), @ (AWOI(IL),IL=1,10), @ (AOE(IL),IL=1,6), @ IPNFI, @ (IPTI(IL),XLI(IL),YLI(IL),IL=1,IPNFI) CALL ROTMAT (AOE(4),AOE(5),AOE(6),SOM1,SPH1,SKP1, @ COM1,CPH1,CKP1,R1) C DX=XJ-AOE(1) DY=YJ-AOE(2) DZ=ZJ-AOE(3) U1=R1(1,1)*DX+R1(1,2)*DY+R1(1,3)*DZ U2=R1(2,1)*DX+R1(2,2)*DY+R1(2,3)*DZ U3=R1(3,1)*DX+R1(3,2)*DY+R1(3,3)*DZ T=U3**2 C C FORM MATRIX (PART.DERIV. WRT OBJECT COORDS) CALL AMAT (AOIUD(3),U1,U2,U3,T,DX,DY,DZ,R1,AA) C EXTRACT PHOTO-COORDS OF PNT LTAB(M,1) ON PHOTO LTAB(M,IP+1) C DO 65 IL=1,IPNFI IF (LTAB(M,1).EQ.IPTI(IL)) THEN XPHOT=XLI(IL) YPHOT=YLI(IL) GO TO 70 ENDIF 65 CONTINUE 70 CONTINUE C FORM MATRIX (MISCLOSURE VECTOR OF PHOTOG OBS OF POINT K) CALL WPMAT (AOIUD,U1,U2,U3,XPHOT,YPHOT,WPA) C = -T CALL TRNSD (AAT,3,AA,2,2,3) C = * PX,PY DO 80 IL=1,3 AAT(IL,1)=AAT(IL,1)*PX AAT(IL,2)=AAT(IL,2)*PY 80 CONTINUE C = * CALL MMULD (C1,3,AAT,3,AA,2,3,2,3) C = * CALL MMULD (AW1,3,AAT,3,WPA,2,3,2,1) C UPDATE AND MATRICES DO 110 JR=1,3 AW(JR)=AW(JR)+AW1(JR) DO 100 IL=1,3 C(IL,JR)=C(IL,JR)+C1(IL,JR) 100 CONTINUE 110 CONTINUE 120 CONTINUE C CONSIDER WEIGHT FACTOR OF THE APPRX. OBJECT COORDS DO 121 IL=1,3 C(IL,IL)=C(IL,IL)+PAUX 121 CONTINUE C C CHECK IF PNT K=LTAB(M,1) IS A CONTROL PNT OR NOT C IF (LTAB(M,12).NE.0) THEN DO 160 IL=1,NGF IF (LTAB(M,1).EQ.IG(IL)) THEN W(1)=WX(IL) W(2)=WY(IL) W(3)=WZ(IL) C FORM (MISCLOSURE VECTOR OF THE CONTROL POINT OBSERVATIONS) IF (LTAB(M,12).EQ.1) THEN WC(1)=XJ-XG(IL) WC(2)=YJ-YG(IL) WC(3)=ZJ-ZG(IL) ELSEIF (LTAB(M,12).EQ.2) THEN WC(1)=XJ-XG(IL) WC(2)=YJ-YG(IL) WC(3)=0.D0 ELSEIF (LTAB(M,12).EQ.3) THEN WC(1)=0.D0 WC(2)=0.D0 WC(3)=ZJ-ZG(IL) ENDIF C COMPLETE AND WITH THE CONTRIBUTION OF THE CONTROL PNT DO 140 JR=1,3 C(JR,JR)=C(JR,JR)+W(JR) AW(JR)=AW(JR)+W(JR)*WC(JR) 140 CONTINUE GO TO 180 ENDIF 160 CONTINUE ENDIF 180 CONTINUE C = -INV CALL SPIN (C,3,3,DET,IDEXP) C C ALL MATRICES FOR PNT K HAVE BEEN COMPUTED. COMPUTE MATRIX SUMMATIONS C ***** FORM ELEMENTS OF MATRIX N11(I,I) AND OF ** C = -T CALL TRNSD (ET,6,E,2,2,6) C = * PX,PY DO 200 IL=1,6 ET(IL,1)=ET(IL,1)*PX ET(IL,2)=ET(IL,2)*PY 200 CONTINUE C = * CALL MMULD (E1,6,ET,6,E,2,6,2,6) C = * CALL MMULD (E2,6,ET,6,A,2,6,2,3) C = * CALL MMULD (E3,6,E2,6,C,3,6,3,3) C = -T CALL TRNSD (E2T,3,E2,6,6,3) C = * CALL MMULD (E4,6,E3,6,E2T,3,6,3,6) C = * CALL MMULD (EW,6,ET,6,WP,2,6,2,1) C = * CALL MMULD (AW2,6,E3,6,AW,3,6,3,1) C UPDATE , AND DO 240 JR=1,6 AW0(JR)=AW0(JR)+EW(JR)-AW2(JR) DO 220 IL=1,6 E0(IL,JR)=E0(IL,JR)+E1(IL,JR)-E4(IL,JR) 220 CONTINUE 240 CONTINUE C CONSIDER IO PARAMETERS IF (ITER.GE.ITRI) THEN C C ************ FORM ELEMENTS OF MATRIX C C =

-TR CALL TRNSD (PT,10,P,2,2,10) C = * PX,PY DO 260 IL=1,10 PT(IL,1)=PT(IL,1)*PX PT(IL,2)=PT(IL,2)*PY 260 CONTINUE C = * CALL MMULD (P1,10,PT,10,E,2,10,2,6) C = * CALL MMULD (P2,10,PT,10,A,2,10,2,3) C = * CALL MMULD (P3,10,P2,10,C,3,10,3,3) C C = -T CALL TRNSD (AT,3,A,2,2,3) C = * PX,PY DO 280 IL=1,3 AT(IL,1)=AT(IL,1)*PX AT(IL,2)=AT(IL,2)*PY 280 CONTINUE C = * CALL MMULD (P4,3,AT,3,E,2,3,2,6) C = * CALL MMULD (P5,10,P3,10,P4,3,10,3,6) C UPDATE WHICH CONTAINS ELEMENTS OF DO 320 JR=1,6 DO 300 IL=1,10 P0(IL,JR)=P0(IL,JR)+P1(IL,JR)-P5(IL,JR) 300 CONTINUE 320 CONTINUE C C C ** FORM ELEMENTS AND OF MATRICES AND ** C C = *

CALL MMULD (Q,10,PT,10,P,2,10,2,10) C = * CALL MMULD (QW,10,PT,10,WP,2,10,2,1) C COUNT THE WEIGHTS FOR THE IO ELEMENTS, Q=Q+WOI & QW=QW+WOI*WI DO 340 IL=1,10 Q(IL,IL)=Q(IL,IL)+WOI(IL) QW(IL)=QW(IL)+WOI(IL)*WI(IL) 340 CONTINUE C = -T CALL TRNSD (Q1,3,P2,10,10,3) C = * CALL MMULD (Q2,10,P3,10,Q1,3,10,3,10) C = * CALL MMULD (QW1,10,P3,10,AW,3,10,3,1) C UPDATE WHICH CONTAIN THE ELEMENTS OF , DO 360 JR=1,10 QW0(JR)=QW0(JR)+QW(JR)-QW1(JR) DO 350 IL=1,10 Q0(IL,JR)=Q0(IL,JR)+Q(IL,JR)-Q2(IL,JR) 350 CONTINUE 360 CONTINUE C FOR IO PARAMETERS ENDIF C 1200 CONTINUE C C END OF PHOTO-PNT CONTRIBUTION OF PHOTO I TO (N11,N21,N22 OF (I,I)) C STORE ,, TO PROPER LOCATIONS IN C SYMMETRIC STORAGE MODE C C IN VIA , IN VIA ME=0 NE=0 I6=6*I II=I6-5 DO 390 IL=II,I6 ME=ME+1 U(IL)=AW0(ME) C UEO(IL)=AW0(ME) DO 380 JR=II,I6 NE=NE+1 C NS(IL,JR)=E0(ME,NE) IF (IL.GE.JR) THEN IPOS=(IL*(IL-1)/2)+JR N(IPOS)=E0(ME,NE) C NEO(IPOS)=E0(ME,NE) ENDIF 380 CONTINUE NE=0 390 CONTINUE C CONSIDER THE IO PARAMETERS IF (ITER.GE.ITRI) THEN C C IN VIA MP=0 NP=0 I6=6*I II=I6-5 J6=6*NPHOT IS=J6+1+10*(I-1) I10=10*I+J6 DO 410 IL=IS,I10 MP=MP+1 DO 400 JR=II,I6 NP=NP+1 C N(IL,JR)=P0(MP,NP) C N(JR,IL)=N(IL,JR) IF (IL.GE.JR) THEN IPOS=(IL*(IL-1)/2)+JR N(IPOS)=P0(MP,NP) ENDIF 400 CONTINUE NP=0 410 CONTINUE C C IN VIA , IN VIA MQ=0 NQ=0 J6=6*NPHOT IS=J6+1+10*(I-1) I10=10*I+J6 DO 430 IL=IS,I10 MQ=MQ+1 U(IL)=QW0(MQ) DO 420 JR=IS,I10 NQ=NQ+1 C N(IL,JR)=Q0(MQ,NQ) IF (IL.GE.JR) THEN IPOS=(IL*(IL-1)/2)+JR N(IPOS)=Q0(MQ,NQ) ENDIF 420 CONTINUE NQ=0 430 CONTINUE C NEW MATRICES NIO C ME=0 C NE=0 C I10=10*I C II=I10-9 C DO 3901 IL=II,I10 C ME=ME+1 C UIO(IL)=QW0(ME) C DO 3801 JR=II,I10 C NE=NE+1 C IF (IL.GE.JR) THEN C IPOS=(IL*(IL-1)/2)+JR C NIO(IPOS)=Q0(ME,NE) C ENDIF C3801 CONTINUE C NE=0 C3901 CONTINUE C END OF IO PARAMETERS ENDIF C ELSEIF (J.LT.I) THEN C C FORM , J C DO 1300 MP=1,JFP IPHOT=0 JPHOT=0 DO 435 IL=1,LTAB(MP,11) C CHECK IF POINT APPEARS ON PHOTOGRAPH I IF (LTAB(MP,IL+1).EQ.I) THEN IPHOT=I GO TO 440 ENDIF 435 CONTINUE 440 CONTINUE DO 450 IL=1,LTAB(MP,11) C CHECK IF POINT APPEARS ON PHOTOGRAPH J IF (LTAB(MP,IL+1).EQ.J) THEN JPHOT=J GO TO 460 ENDIF 450 CONTINUE 460 CONTINUE C IF (IPHOT.NE.0.AND.JPHOT.NE.0) THEN C POINT APPEARS ON BOTH PHOTOS I AND J C GET PHOTO-RECORD I READ (UNIT=IEX1,REC=IPHOT,ERR=9900) IPH, @ JPH, @ (OIUD(K),K=1,10), @ (OI(K),K=1,10), @ (WOI(K),K=1,10), @ (OE(K),K=1,6), @ IPNF, @ (IPT(K),XL(K),YL(K),K=1,IPNF) CALL ROTMAT (OE(4),OE(5),OE(6),SOM,SPH,SKP,COM,CPH,CKP,R) C C EXTRACT PHOTO-COORDS OF PNT LTAB(MP,1) ON PHOTO I=IPHOT C DO 465 IL=1,IPNF IF (LTAB(MP,1).EQ.IPT(IL)) THEN XPHOT=XL(IL) YPHOT=YL(IL) GO TO 467 ENDIF 465 CONTINUE 467 CONTINUE C GET APPRX. OBJECT COORDS OF PNT LTAB(MP,1) DO 470 L=1,JFP READ(UNIT=IEX2,REC=L,ERR=3500) IDO,XJ,YJ,ZJ IF (LTAB(MP,1).EQ.IDO) THEN DX=XJ-OE(1) DY=YJ-OE(2) DZ=ZJ-OE(3) GO TO 480 ENDIF 470 CONTINUE 480 CONTINUE U1=R(1,1)*DX+R(1,2)*DY+R(1,3)*DZ U2=R(2,1)*DX+R(2,2)*DY+R(2,3)*DZ U3=R(3,1)*DX+R(3,2)*DY+R(3,3)*DZ T=U3**2 C FORM MATRIX (PART.DERIV. WRT EO ELEMENTS) CALL EMAT (OIUD(3),U1,U2,U3,T,DX,DY,DZ,R,SOM,SPH,SKP, @ COM,CPH,CKP,E) C C FORM MATRIX (PART.DERIV. WRT OBJECT COORDS) CALL AMAT (OIUD(3),U1,U2,U3,T,DX,DY,DZ,R,A) C C FORM

MATRIX (PART.DERIV. WRT IO ELEMENTS) CALL PMAT (OIUD,U1,U2,U3,XPHOT,YPHOT,P) C C FORM MATRIX FOR THE LTAB(MP,1) PNT. USE MATRIX INSTEAD OF C INITIALIAZE FOR THE LTAB(MP,1) PNT DO 500 JR=1,3 DO 490 IL=1,3 C(IL,JR)=0.D0 490 CONTINUE 500 CONTINUE C GET THE PHOTOS ON WHICH PNT LTAB(MP,1) APPEARS DO 540 IP=1,LTAB(MP,11) READ (UNIT=IEX1,REC=LTAB(MP,IP+1),ERR=9900) IPHI, @ JPHI, @ (AOIUD(IL),IL=1,10), @ (AOI(IL),IL=1,10), @ (AWOI(IL),IL=1,10), @ (AOE(IL),IL=1,6) CALL ROTMAT (AOE(4),AOE(5),AOE(6),SOM1,SPH1,SKP1, @ COM1,CPH1,CKP1,R1) C DX=XJ-AOE(1) DY=YJ-AOE(2) DZ=ZJ-AOE(3) U1=R1(1,1)*DX+R1(1,2)*DY+R1(1,3)*DZ U2=R1(2,1)*DX+R1(2,2)*DY+R1(2,3)*DZ U3=R1(3,1)*DX+R1(3,2)*DY+R1(3,3)*DZ T=U3**2 C C FORM MATRIX (PART.DERIV. WRT OBJECT COORDS) CALL AMAT (AOIUD(3),U1,U2,U3,T,DX,DY,DZ,R1,AA) C C = -T CALL TRNSD (AAT,3,AA,2,2,3) C = * PX,PY DO 510 IL=1,3 AAT(IL,1)=AAT(IL,1)*PX AAT(IL,2)=AAT(IL,2)*PY 510 CONTINUE C = * CALL MMULD (C1,3,AAT,3,AA,2,3,2,3) C UPDATE MATRIX DO 530 JR=1,3 DO 520 IL=1,3 C(IL,JR)=C(IL,JR)+C1(IL,JR) 520 CONTINUE 530 CONTINUE 540 CONTINUE DO 541 IL=1,3 C(IL,IL)=C(IL,IL)+PAUX 541 CONTINUE C C CHECK IF PNT LTAB(MP,1) IS A CONTROL PNT OR NOT C IF (LTAB(MP,12).NE.0) THEN DO 560 IL=1,NGF IF (LTAB(MP,1).EQ.IG(IL)) THEN W(1)=WX(IL) W(2)=WY(IL) W(3)=WZ(IL) C COMPLETE WITH THE CONTRIBUTION OF THE CONTROL PNT DO 550 JR=1,3 C(JR,JR)=C(JR,JR)+W(JR) 550 CONTINUE GO TO 580 ENDIF 560 CONTINUE ENDIF 580 CONTINUE C = -INV CALL SPIN (C,3,3,DET,IDEXP) C = -TR FORM MATRICES FOR N11,N21,N22 CALL TRNSD (ET,6,E,2,2,6) C = * PX,PY DO 590 IL=1,6 ET(IL,1)=ET(IL,1)*PX ET(IL,2)=ET(IL,2)*PY 590 CONTINUE C = * CALL MMULD (E2,6,ET,6,A,2,6,2,3) C = * ******** FOR N11 ********* CALL MMULD (E3,6,E2,6,C,3,6,3,3) C =

-TR CALL TRNSD (PT,10,P,2,2,10) C = * PX,PY DO 610 IL=1,10 PT(IL,1)=PT(IL,1)*PX PT(IL,2)=PT(IL,2)*PY 610 CONTINUE C = * CALL MMULD (P2,10,PT,10,A,2,10,2,3) C = * ****** FOR N21 & N22 ******** CALL MMULD (P3,10,P2,10,C,3,10,3,3) C END OF MATRIX FORMATIONS FOR PNT LTAB(MP,1) ON PHOTO I C GET PHOTO-RECORD J READ (UNIT=IEX1,REC=JPHOT,ERR=9900) IPH, @ JPH, @ (OIUD(K),K=1,10), @ (OI(K),K=1,10), @ (WOI(K),K=1,10), @ (OE(K),K=1,6), @ IPNF, @ (IPT(K),XL(K),YL(K),K=1,IPNF) CALL ROTMAT (OE(4),OE(5),OE(6),SOM,SPH,SKP,COM,CPH,CKP,R) C C EXTRACT PHOTO-COORDS OF PNT LTAB(MP,1) ON PHOTO J=JPHOT C DO 615 IL=1,IPNF IF (LTAB(MP,1).EQ.IPT(IL)) THEN XPHOT=XL(IL) YPHOT=YL(IL) GO TO 617 ENDIF 615 CONTINUE 617 CONTINUE DX=XJ-OE(1) DY=YJ-OE(2) DZ=ZJ-OE(3) U1=R(1,1)*DX+R(1,2)*DY+R(1,3)*DZ U2=R(2,1)*DX+R(2,2)*DY+R(2,3)*DZ U3=R(3,1)*DX+R(3,2)*DY+R(3,3)*DZ T=U3**2 C FORM MATRIX (PART.DERIV. WRT EO ELEMENTS) CALL EMAT (OIUD(3),U1,U2,U3,T,DX,DY,DZ,R,SOM,SPH,SKP, @ COM,CPH,CKP,E) C C FORM MATRIX (PART.DERIV. WRT OBJECT COORDS) CALL AMAT (OIUD(3),U1,U2,U3,T,DX,DY,DZ,R,A) C C FORM

MATRIX (PART.DERIV. WRT IO ELEMENTS) CALL PMAT (OIUD,U1,U2,U3,XPHOT,YPHOT,P) C = -TR CALL TRNSD (AT,3,A,2,2,3) C = * PX,PY DO 620 IL=1,3 AT(IL,1)=AT(IL,1)*PX AT(IL,2)=AT(IL,2)*PY 620 CONTINUE C = * ******** FOR N11,N21 ****** CALL MMULD (P4,3,AT,3,E,2,3,2,6) C = *

******** FOR N22 ******** CALL MMULD (Q1,3,AT,3,P,2,3,2,10) C ALL NECESSARY MATRIX COMPONENTS ARE FORMED. FORM N1, J = * CALL MMULD (E4,6,E3,6,P4,3,6,3,6) C UPDATE WHICH CONTAINS THE ELEMENTS OF , J = * FORM N21, J WHICH CONTAINS THE ELEMENTS OF , J = * FORM N22, J WHICH CONTAINS THE ELEMENTS OF , J,J VIA ME=0 NE=0 JE=6*J JB=JE-5 IE=6*I IB=IE-5 DO 700 IL=IB,IE ME=ME+1 DO 690 JR=JB,JE NE=NE+1 C NS(IL,JR)=E0(ME,NE) C NS(JR,IL)=NS(IL,JR) IF (IL.GE.JR) THEN IPOS=(IL*(IL-1)/2)+JR N(IPOS)=E0(ME,NE) C NEO(IPOS)=E0(ME,NE) ENDIF 690 CONTINUE NE=0 700 CONTINUE C CONSIDER IO PARAMETERS IF (ITER.GE.ITRI) THEN C C ,J VIA MP=0 NP=0 JE=6*J JB=JE-5 I6=6*NPHOT IB=I6+1+10*(I-1) IE=10*I+I6 DO 720 IL=IB,IE MP=MP+1 DO 710 JR=JB,JE NP=NP+1 C N(IL,JR)=P0(MP,NP) C N(JR,IL)=N(IL,JR) IF (IL.GE.JR) THEN IPOS=(IL*(IL-1)/2)+JR N(IPOS)=P0(MP,NP) ENDIF 710 CONTINUE NP=0 720 CONTINUE C C ,J VIA MQ=0 NQ=0 I6=6*NPHOT JB=I6+1+10*(J-1) JE=10*J+I6 IB=I6+1+10*(I-1) IE=10*I+I6 DO 740 IL=IB,IE MQ=MQ+1 DO 730 JR=JB,JE NQ=NQ+1 C N(IL,JR)=Q0(MQ,NQ) C N(JR,IL)=N(IL,JR) IF (IL.GE.JR) THEN IPOS=(IL*(IL-1)/2)+JR N(IPOS)=Q0(MQ,NQ) ENDIF 730 CONTINUE NQ=0 740 CONTINUE C C ME=0 C NE=0 C JE=10*J C JB=JE-9 C IE=10*I C IB=IE-9 C DO 7001 IL=IB,IE C ME=ME+1 C DO 6901 JR=JB,JE C NE=NE+1 C N(IL,JR)=Q0(ME,NE) C N(JR,IL)=N(IL,JR) C IF (IL.GE.JR) THEN C IPOS=(IL*(IL-1)/2)+JR C NIO(IPOS)=Q0(ME,NE) C ENDIF C6901 CONTINUE C NE=0 C7001 CONTINUE C END OF IO PARAMETERS ENDIF C ELSEIF (J.GT.I) THEN C CONSIDER IO PARAMETERS IF (ITER.GE.ITRI) THEN C C FORM , J>I (ABOVE DIAGONAL ELEMENTS) C FOR ALL POINTS IN C DO 1500 MP=1,JFP IPHOT=0 JPHOT=0 DO 750 IL=1,LTAB(MP,11) C CHECK IF POINT APPEARS ON PHOTOGRAPH I IF (LTAB(MP,IL+1).EQ.I) THEN IPHOT=I GO TO 760 ENDIF 750 CONTINUE 760 CONTINUE DO 770 IL=1,LTAB(MP,11) C CHECK IF POINT APPEARS ON PHOTOGRAPH J IF (LTAB(MP,IL+1).EQ.J) THEN JPHOT=J GO TO 780 ENDIF 770 CONTINUE 780 CONTINUE C IF (IPHOT.NE.0.AND.JPHOT.NE.0) THEN C POINT APPEARS ON BOTH PHOTOS I AND J C GET PHOTO-RECORD I READ (UNIT=IEX1,REC=IPHOT,ERR=9900) IPH, @ JPH, @ (OIUD(K),K=1,10), @ (OI(K),K=1,10), @ (WOI(K),K=1,10), @ (OE(K),K=1,6), @ IPNF, @ (IPT(K),XL(K),YL(K),K=1,IPNF) CALL ROTMAT (OE(4),OE(5),OE(6),SOM,SPH,SKP,COM,CPH,CKP,R) C C EXTRACT PHOTO-COORDS OF PNT LTAB(MP,1) ON PHOTO I=IPHOT C DO 785 IL=1,IPNF IF (LTAB(MP,1).EQ.IPT(IL)) THEN XPHOT=XL(IL) YPHOT=YL(IL) GO TO 787 ENDIF 785 CONTINUE 787 CONTINUE C DO 790 L=1,JFP READ(UNIT=IEX2,REC=L,ERR=3500) IDO,XJ,YJ,ZJ IF (LTAB(MP,1).EQ.IDO) THEN DX=XJ-OE(1) DY=YJ-OE(2) DZ=ZJ-OE(3) GO TO 800 ENDIF 790 CONTINUE 800 CONTINUE U1=R(1,1)*DX+R(1,2)*DY+R(1,3)*DZ U2=R(2,1)*DX+R(2,2)*DY+R(2,3)*DZ U3=R(3,1)*DX+R(3,2)*DY+R(3,3)*DZ T=U3**2 C FORM MATRIX (PART.DERIV. WRT OBJECT COORDS) CALL AMAT (OIUD(3),U1,U2,U3,T,DX,DY,DZ,R,A) C C FORM

MATRIX (PART.DERIV. WRT IO ELEMENTS) CALL PMAT (OIUD,U1,U2,U3,XPHOT,YPHOT,P) C C FORM MATRIX FOR THE LTAB(MP,1) PNT. USE MATRIX INSTEAD OF C INITIALIAZE FOR THE LTAB(MP,1) PNT DO 820 JR=1,3 DO 810 IL=1,3 C(IL,JR)=0.D0 810 CONTINUE 820 CONTINUE C GET THE PHOTOS ON WHICH PNT LTAB(MP,1) APPEARS DO 860 IP=1,LTAB(MP,11) READ (UNIT=IEX1,REC=LTAB(MP,IP+1),ERR=9900) IPHI, @ JPHI, @ (AOIUD(IL),IL=1,10), @ (AOI(IL),IL=1,10), @ (AWOI(IL),IL=1,10), @ (AOE(IL),IL=1,6) CALL ROTMAT (AOE(4),AOE(5),AOE(6),SOM1,SPH1,SKP1, @ COM1,CPH1,CKP1,R1) C DX=XJ-AOE(1) DY=YJ-AOE(2) DZ=ZJ-AOE(3) U1=R1(1,1)*DX+R1(1,2)*DY+R1(1,3)*DZ U2=R1(2,1)*DX+R1(2,2)*DY+R1(2,3)*DZ U3=R1(3,1)*DX+R1(3,2)*DY+R1(3,3)*DZ T=U3**2 C C FORM MATRIX (PART.DERIV. WRT OBJECT COORDS) CALL AMAT (AOIUD(3),U1,U2,U3,T,DX,DY,DZ,R1,AA) C C = -T CALL TRNSD (AAT,3,AA,2,2,3) C = * PX,PY DO 830 IL=1,3 AAT(IL,1)=AAT(IL,1)*PX AAT(IL,2)=AAT(IL,2)*PY 830 CONTINUE C = * CALL MMULD (C1,3,AAT,3,AA,2,3,2,3) C UPDATE MATRIX DO 850 JR=1,3 DO 840 IL=1,3 C(IL,JR)=C(IL,JR)+C1(IL,JR) 840 CONTINUE 850 CONTINUE 860 CONTINUE DO 861 IL=1,3 C(IL,IL)=C(IL,IL)+PAUX 861 CONTINUE C C CHECK IF PNT LTAB(MP,1) IS A CONTROL PNT OR NOT C IF (LTAB(MP,12).NE.0) THEN DO 880 IL=1,NGF IF (LTAB(MP,1).EQ.IG(IL)) THEN W(1)=WX(IL) W(2)=WY(IL) W(3)=WZ(IL) C COMPLETE WITH THE CONTRIBUTION OF THE CONTROL PNT DO 870 JR=1,3 C(JR,JR)=C(JR,JR)+W(JR) 870 CONTINUE GO TO 890 ENDIF 880 CONTINUE ENDIF 890 CONTINUE C = -INV CALL SPIN (C,3,3,DET,IDEXP) C FORM MATRICES FOR N21 C =

-TR CALL TRNSD (PT,10,P,2,2,10) C = * PX,PY DO 900 IL=1,10 PT(IL,1)=PT(IL,1)*PX PT(IL,2)=PT(IL,2)*PY 900 CONTINUE C = * CALL MMULD (P2,10,PT,10,A,2,10,2,3) C = * ****** FOR N21 ******** CALL MMULD (P3,10,P2,10,C,3,10,3,3) C END OF MATRIX FORMATIONS FOR PNT LTAB(MP,1) ON PHOTO I C GET PHOTO-RECORD J READ (UNIT=IEX1,REC=JPHOT,ERR=9900) IPH, @ JPH, @ (OIUD(K),K=1,10), @ (OI(K),K=1,10), @ (WOI(K),K=1,10), @ (OE(K),K=1,6), @ IPNF, @ (IPT(K),XL(K),YL(K),K=1,IPNF) CALL ROTMAT (OE(4),OE(5),OE(6),SOM,SPH,SKP,COM,CPH,CKP,R) C DX=XJ-OE(1) DY=YJ-OE(2) DZ=ZJ-OE(3) U1=R(1,1)*DX+R(1,2)*DY+R(1,3)*DZ U2=R(2,1)*DX+R(2,2)*DY+R(2,3)*DZ U3=R(3,1)*DX+R(3,2)*DY+R(3,3)*DZ T=U3**2 C FORM MATRIX (PART.DERIV. WRT EO ELEMENTS) CALL EMAT (OIUD(3),U1,U2,U3,T,DX,DY,DZ,R,SOM,SPH,SKP, @ COM,CPH,CKP,E) C C FORM MATRIX (PART.DERIV. WRT OBJECT COORDS) CALL AMAT (OIUD(3),U1,U2,U3,T,DX,DY,DZ,R,A) C C = -TR CALL TRNSD (AT,3,A,2,2,3) C = * PX,PY DO 910 IL=1,3 AT(IL,1)=AT(IL,1)*PX AT(IL,2)=AT(IL,2)*PY 910 CONTINUE C = * ******** FOR N21 ****** CALL MMULD (P4,3,AT,3,E,2,3,2,6) C ALL NECESSARY MATRIX COMPONENTS ARE FORMED C = * FORM N21, J>I CALL MMULD (P5,10,P3,10,P4,3,10,3,6) C UPDATE WHICH CONTAINS THE ELEMENTS OF , J>I DO 930 JR=1,6 DO 920 IL=1,10 P0(IL,JR)=P0(IL,JR)-P5(IL,JR) 920 CONTINUE 930 CONTINUE ENDIF 1500 CONTINUE C ,J>I IN VIA MP=0 NP=0 JE=6*J JB=JE-5 I6=6*NPHOT IB=I6+1+10*(I-1) IE=10*I+I6 DO 960 IL=IB,IE MP=MP+1 DO 950 JR=JB,JE NP=NP+1 C N(IL,JR)=P0(MP,NP) C N(JR,IL)=N(IL,JR) IF (IL.GE.JR) THEN IPOS=(IL*(IL-1)/2)+JR N(IPOS)=P0(MP,NP) ENDIF 950 CONTINUE NP=0 960 CONTINUE C END OF IO PARAMETERS ENDIF C ENDIF 2000 CONTINUE 3000 CONTINUE C CALCULATION OF ACTUAL DIMENSIONS MIO=10*NPHOT MEO=6*NPHOT MPR=MEO+MIO C IDMNE=(MEO*(MEO+1))/2 C IDMNI=(MIO*(MIO+1))/2 C WRITE(IW,*) ' ELEMENTS OF N' C DO 8888 IL=1,MPR C DO 7777 JR=1,IL C IPOS=(IL*(IL-1)/2)+JR C WRITE(IW,*) IL,JR,N(IPOS) C7777 CONTINUE C8888 CONTINUE C C WRITE(6,*) ' NS MATRIX ******' C CALL MOUTD (NS,12,12,12) C WRITE(6,*) ' UEO MATRIX *****' C CALL MOUTD (UEO,30,12,1) C CALL SPIN (NS,12,12,DET,IDEXP) C WRITE(6,*) ' DET= ',DET,' IDEXP= ',IDEXP C CALL MMULD (XX,12,NS,12,UEO,30,12,12,1) C WRITE(6,*) ' XX MATRIX *****' C CALL MOUTD (XX,12,12,1) C RETURN C C WRITE(6,*) ' CORRECTIONS FOR EO & IO' C CALL SYMSYS (N,1,MPR,U,JDMN1,IDGT,D1,D2,IERR) C CALL MOUTD(U,JDMN1,MPR,1) C IF (ITER.EQ.10) THEN C GET THE SOLUTION (CORRECTIONS), IN , FOR THE CAMERA EO WRITE(IW,*) ' CORRECTIONS FOR EO' C C CALL LEQT2P (NEO,1,MEO,JDMN1,UEO,IDGT,D1,D2,WKAREA,IERR) C CALL SYMSYS (NEO,1,MEO,UEO,JDMN1,IDGT,D1,D2,IERR) C IF (ITER.GE.ITRI) THEN JU=MPR ELSE JU=MEO ENDIF C CALL SYMSYS (N,1,JU,U,NUDIM,IDGT,D1,D2,IERR) C WRITE(IW,*) ' D1= ',D1 WRITE(IW,*) ' D2= ',D2 WRITE(IW,*) ' IER= ',IERR IF (IERR.NE.0) RETURN C CALL MOUTD (UEO,30,MEO,1) CALL MOUTD (U,NUDIM,JU,1) C GET SOLUTION (CORRECTIONS), IN , FOR THE CAMERA IO C C IF (ITER.GE.ITRI) THEN C CALL IOCOR (IDMNI,MIO,NIO,UIO,IERR) C ENDIF C IF (IERR.NE.0) RETURN C C ENDIF C C GET SOLUTION (CORRECTIONS) FOR THE APPRX. OBJECT COORDINATES & UPDATE C C CALL OBCOR1(MEO,MIO,WP,AW,AW1,W,WC,A,AT,C,C1,E,P4,P,Q1,UEO,UIO) C CALL OBCOR (JU,WP,AW,AW1,W,WC,A,AT,C,C1,E,P4,P,Q1,U,QVB) C C UPDATE EO/IO PARAMETERS FOR EACH PHOTOGRAPH WRITE(IW,*) ' IO/EO PHOTO-PARAMETERS' C DO 990 IL=1,NPHOT READ (UNIT=IEX1,REC=IL,ERR=9900) IPH, @ JPH, @ (OIUD(JR),JR=1,10), @ (OI(JR),JR=1,10), @ (WOI(JR),JR=1,10), @ (OE(JR),JR=1,6), @ IPNF, @ (IPT(JR),XL(JR),YL(JR),JR=1,IPNF) C UPDATE EO PARAMETERS. ADD THE CORRESPONDING CORRECTIONS FROM JCE=(IL-1)*6 DO 980 JR=1,6 JCE=JCE+1 OE(JR)=OE(JR)-U(JCE) C OE(JR)=OE(JR)-UEO(JCE) 980 CONTINUE IF (ITER.GE.ITRI) THEN C UPDATE IO PARAMETERS. ADD THE CORRESPONDING CORRECTIONS FROM JCI=(IL-1)*10+NPHOT*6 C JCI=(IL-1)*10 DO 985 JR=1,10 JCI=JCI+1 OIUD(JR)=OIUD(JR)-U(JCI) C OIUD(JR)=OIUD(JR)-UIO(JCI) 985 CONTINUE ENDIF C STORE THE UPDATED EO/IO PARAMETERS OF EACH PHOTOGRAPH ON WRITE (UNIT=IEX1,REC=IL,ERR=9900) IPH, @ JPH, @ (OIUD(JR),JR=1,10), @ (OI(JR),JR=1,10), @ (WOI(JR),JR=1,10), @ (OE(JR),JR=1,6), @ IPNF, @ (IPT(JR),XL(JR),YL(JR),JR=1,IPNF) WRITE(IW,*) JPH,IPH WRITE(IW,*) (OIUD(JR),JR=1,10) WRITE(IW,*) (OE(JR),JR=1,6) 990 CONTINUE C C UPDATE NUMBER OF ITERATIONS ITER=ITER+1 C C COMPARE ITER WITH ITRF IF (ITER.LE.ITRF) GO TO 6010 C COMPUTE PHOTO-RESIDUALS CALL PHORES (W1,QVP,QVI) C COMPUTE RESIDUALS AT CHECK POINTS CALL CEKRES C COMPUTE RESIDUALS AT CONTROL POINTS CALL CEKCON (W,WC,QVP,QVI,QVB,SOC) C COMPUTE VARIANCE-COVARIANCE MATRIX OF OBJECT POINTS CALL COVAR (JU,NDIM,N,SOC,A,AT,C,C1,E,P4,P,Q1,E2T,ET,E2,P1,Q2, @ PT,P2) RETURN 3500 WRITE(IW,3501) 3501 FORMAT(' ',10X,'ERROR RELATED TO FILE ',/) RETURN 9900 WRITE(IW,9901) 9901 FORMAT(' ',10X,'ERROR RELATED TO FILE ',/) RETURN END C C ******************************************************************* SUBROUTINE EMAT (F,U1,U2,U3,T,DX,DY,DZ,R,SOM,SPH,SKP,COM,CPH, @ CKP,E) C C FORMS MATRIX PARTIAL DERIVATIVES WRT EO ELEMENTS C * WRITTEN BY: COSTAS ARMENAKIS ON JULY 1986 * C ******************************************************************* C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION R(3,3),E(2,6) C T1=F/T T2=(CPH*DX+R(3,1)*SOM*DY-R(3,1)*COM*DZ)*U1 T3=R(3,1)*U1 T4=R(3,2)*U1 T5=R(3,3)*U1 C E(1,1)=T1*(R(1,1)*U3-T3) E(1,2)=T1*(R(1,2)*U3-T4) E(1,3)=T1*(R(1,3)*U3-T5) E(1,4)=((R(1,3)*DY-R(1,2)*DZ)*U3+(-R(3,3)*DY+R(3,2)*DZ)*U1)*T1 E(1,5)=((R(3,1)*CKP*DX-R(1,1)*SOM*DY+R(1,1)*COM*DZ)*U3+T2)*T1 E(1,6)=(-R(2,1)*DX-R(2,2)*DY-R(2,3)*DZ)*U3*T1 C E(2,1)=T1*(R(2,1)*U3-T3) E(2,2)=T1*(R(2,2)*U3-T4) E(2,3)=T1*(R(2,3)*U3-T5) E(2,4)=((R(2,3)*DY-R(2,2)*DZ)*U3+(-R(3,3)*DY+R(3,2)*DZ)*U1)*T1 E(2,5)=((-R(3,1)*SKP*DX-R(2,1)*SOM*DY+R(2,1)*COM*DZ)*U3+T2)*T1 E(2,6)=(R(1,1)*DX+R(1,2)*DY+R(1,3)*DZ)*U3*T1 RETURN END C C ******************************************************************* SUBROUTINE AMAT (F,U1,U2,U3,T,DX,DY,DZ,R,A) C C FORMS MATRIX PARTIAL DERIVATIVES WRT OBJECT COORDS C * WRITTEN BY: COSTAS ARMENAKIS ON JULY 1986 * C ******************************************************************* C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION R(3,3),A(2,3) C T0=F/T T1=R(3,1)*U1 T2=R(3,2)*U1 T3=R(3,3)*U1 C A(1,1)=T0*(-R(1,1)*U3+T1) A(1,2)=T0*(-R(1,2)*U3+T2) A(1,3)=T0*(-R(1,3)*U3+T3) C A(2,1)=T0*(-R(2,1)*U3+T1) A(2,2)=T0*(-R(2,2)*U3+T2) A(2,3)=T0*(-R(2,3)*U3+T3) RETURN END C C ******************************************************************* SUBROUTINE PMAT (O,U1,U2,U3,X,Y,P) C C FORMS MATRIX

PARTIAL DERIVATIVES WRT IO ELEMENTS C * WRITTEN BY: COSTAS ARMENAKIS ON JULY 1986 * C ******************************************************************* C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION O(10),P(2,10) C SX=X-O(1) SY=Y-O(2) SX2=SX**2 SY2=SY**2 S2=SX2+SY2 S4=S2**2 S6=S2**3 B1=O(6)+2.D0*O(7)*S2+3.D0*O(8)*S4 B2=O(6)+O(7)*S2+O(8)*S4 C P(1,1)=1.D0-O(4)-2.D0*SX2*B1-S2*B2-2.D0*O(9)*SY-6.D0*O(10)*SX P(1,2)=-O(5)-2.D0*SX*SY*B1-2.D0*O(9)*SX-2.D0*O(10)*SY P(1,3)=-U1/U3 P(1,4)=SX P(1,5)=SY P(1,6)=S2*SX P(1,7)=S4*SX P(1,8)=S6*SX P(1,9)=2.D0*SX*SY P(1,10)=S2+2.D0*SX2 C P(2,1)=P(1,2) P(2,2)=1.D0+O(4)-2.D0*SY2*B1-S2*B2-6.D0*O(9)*SY-2.D0*O(10)*SX P(2,3)=-U2/U3 P(2,4)=-P(1,5) P(2,5)=P(1,4) P(2,6)=S2*SY P(2,7)=S4*SY P(2,8)=S6*SY P(2,9)=S2+2.D0*SY2 P(2,10)=P(1,9) RETURN END C C ******************************************************************* SUBROUTINE WPMAT (O,U1,U2,U3,X,Y,WP) C C FORMS MISCLOSURE VECTOR OF THE PHOTOGRAMMETRIC OBSERVATIONS C * WRITTEN BY: COSTAS ARMENAKIS ON JULY 1986 * C ******************************************************************* C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION O(10),WP(2) C SX=X-O(1) SY=Y-O(2) SXY=SX*SY SX2=SX**2 SY2=SY**2 S2=SX2+SY2 S4=S2**2 S6=S2**3 B=O(6)*S2+O(7)*S4+O(8)*S6 C DX=O(4)*SX+O(5)*SY+B*SX+2.D0*O(9)*SXY+O(10)*(S2+2.D0*SX2) DY=-O(4)*SY+O(5)*SX+B*SY+O(9)*(S2+2.D0*SY2)+2.D0*O(10)*SXY T=O(3)/U3 WP(1)=O(1)+DX-T*U1-X WP(2)=O(2)+DY-T*U2-Y RETURN END C C *********************************************************************** SUBROUTINE WIMAT (OIUD,OI,WI) C C FORMS MISCLOSURE VECTOR OF THE INTERIOR ORIENTATION OBSERVATIONS C * WRITTEN BY: COSTAS ARMENAKIS ON JULY 1986 * C *********************************************************************** C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION OIUD(10),OI(10),WI(10) C DO 20 IL=1,10 WI(IL)=OIUD(IL)-OI(IL) 20 CONTINUE RETURN END C C *********************************************************************** C SUBROUTINE OBCOR1(MEO,MIO,WP,AW,AW1,W,WC,A,AT,C,C1,E,P4,P,Q1,UEO, C @ UIO) SUBROUTINE OBCOR (MPR,WP,AW,AW1,W,WC,A,AT,C,C1,E,P4,P,Q1,U,QVB) C C COMPUTES THE CORRECTIONS (SOLUTION) TO OBJECT COORDINATES C * WRITTEN BY: COSTAS ARMENAKIS ON AUGUST 1986 * C *********************************************************************** C IMPLICIT REAL*8 (A-H,O-Z) REAL*8 WP(2),AW(3),P4E0(3),Q1I0(3),AW1(3),W(3),WC(3),P4E(3), @ Q1I(3),A(2,3),AT(3,2),C(3,3),C1(3,3),E(2,6),P4(3,6), @ P(2,10),Q1(3,10),U(MPR) C @ P(2,10),Q1(3,10),UEO(MEO),UIO(MIO) C COMMON /ARR/ XL(50),YL(50),DID(100),XG(100),YG(100),ZG(100), @ WX(100),WY(100),WZ(100),XK(100),YK(100),ZK(100), @ IPT(50),IG(100),IK(100), @ INDEX(100),IAPN(450),IIPH(450),LTAB(1000,12) COMMON /CAM/ OE(6),OI(10),OIUD(10),WOI(10),R(3,3) COMMON /DEV/ IR,IW,IEX1,IEX2,IBYT1,IBYT2 COMMON /VAR/ PX,PY,PAUX,NPHOT,JFP,NGF,ITRF,ITRI,ITER,NC,LDEF C QVB=QUADRATIC VALUE QVB = 0.D0 WRITE(IW,*) ' OBJECT COORDS AND CORRECTIONS' C FOR EACH POINT M IN C DO 1000 M=1,JFP C C INITIALIZE MATRICES TO ZERO C DO 20 IL=1,3 AW(IL)=0.D0 P4E0(IL)=0.D0 Q1I(IL)=0.D0 Q1I0(IL)=0.D0 DO 10 JR=1,3 C(JR,IL)=0.D0 10 CONTINUE 20 CONTINUE C C GET THE APPRX. COORDS OF PNT M C READ(UNIT=IEX2,REC=M,ERR=5000) IDO,XJ,YJ,ZJ C C GET PHOTOS ON WHICH PNT M APPEARS C DO 190 IP=1,LTAB(M,11) READ (UNIT=IEX1,REC=LTAB(M,IP+1),ERR=6000) IPH, @ JPH, @ (OIUD(K),K=1,10), @ (OI(K),K=1,10), @ (WOI(K),K=1,10), @ (OE(K),K=1,6), @ IPNF, @ (IPT(K),XL(K),YL(K),K=1,IPNF) CALL ROTMAT (OE(4),OE(5),OE(6),SOM,SPH,SKP,COM,CPH,CKP,R) C DX=XJ-OE(1) DY=YJ-OE(2) DZ=ZJ-OE(3) U1=R(1,1)*DX+R(1,2)*DY+R(1,3)*DZ U2=R(2,1)*DX+R(2,2)*DY+R(2,3)*DZ U3=R(3,1)*DX+R(3,2)*DY+R(3,3)*DZ T=U3**2 C EXTRACT PHOTO-COORDS OF PNT LTAB(M,1) ON PHOTO LTAB(M,IP+1) C DO 50 IL=1,IPNF IF (LTAB(M,1).EQ.IPT(IL)) THEN XPHOT=XL(IL) YPHOT=YL(IL) GO TO 60 ENDIF 50 CONTINUE 60 CONTINUE C FORM MATRIX (PART.DERIV. WRT EO ELEMENTS) CALL EMAT (OIUD(3),U1,U2,U3,T,DX,DY,DZ,R,SOM,SPH,SKP, @ COM,CPH,CKP,E) C C FORM MATRIX (PART.DERIV. WRT OBJECT COORDS) CALL AMAT (OIUD(3),U1,U2,U3,T,DX,DY,DZ,R,A) C C FORM

MATRIX (PART.DERIV. WRT IO ELEMENTS) CALL PMAT (OIUD,U1,U2,U3,XPHOT,YPHOT,P) C FORM MATRIX (MISCLOSURE VECTOR OF PHOTOGRAMMETRIC OBSERVATIONS) CALL WPMAT (OIUD,U1,U2,U3,XPHOT,YPHOT,WP) C C = - TRANSP CALL TRNSD (AT,3,A,2,2,3) C = * PX,PY DO 80 IL=1,3 AT(IL,1)=AT(IL,1)*PX AT(IL,2)=AT(IL,2)*PY 80 CONTINUE C = * CALL MMULD (C1,3,AT,3,A,2,3,2,3) C = * CALL MMULD (P4,3,AT,3,E,2,3,2,6) C = * OF PHOTO LTAB(M,IP+1). CONSIDER EO OF PHOTO DO 120 IL=1,3 JCE=(LTAB(M,IP+1)-1)*6 P4E(IL)=0.D0 DO 100 JR=1,6 JCE=JCE+1 P4E(IL)=P4E(IL)+P4(IL,JR)*(-U(JCE)) C P4E(IL)=P4E(IL)+P4(IL,JR)*(-UEO(JCE)) 100 CONTINUE 120 CONTINUE C IF (ITER.GE.ITRI) THEN C = *

CALL MMULD (Q1,3,AT,3,P,2,3,2,10) C = * OF PHOTO LTAB(M,IP+1). CONSIDER IO OF PHOTO DO 150 IL=1,3 C JCI=(LTAB(M,IP+1)-1)*10 JCI=(LTAB(M,IP+1)-1)*10+NPHOT*6 Q1I(IL)=0.D0 DO 140 JR=1,10 JCI=JCI+1 Q1I(IL)=Q1I(IL)+Q1(IL,JR)*(-U(JCI)) C Q1I(IL)=Q1I(IL)+Q1(IL,JR)*(-UIO(JCI)) 140 CONTINUE 150 CONTINUE ENDIF C = * CALL MMULD (AW1,3,AT,3,WP,2,3,2,1) C UPDATE MATRICES FOR EACH PHOTO DO 180 IL=1,3 AW(IL)=AW(IL)+AW1(IL) P4E0(IL)=P4E0(IL)+P4E(IL) Q1I0(IL)=Q1I0(IL)+Q1I(IL) DO 160 JR=1,3 C(JR,IL)=C(JR,IL)+C1(JR,IL) 160 CONTINUE 180 CONTINUE C 190 CONTINUE C CONSIDER THE WEIGHT FACTOR OF THE APPRX. OBJECT COORDS DO 200 IL=1,3 C(IL,IL)=C(IL,IL)+PAUX 200 CONTINUE C C CHECK IF PNT M=LTAB(M,1) IS A CONTROL PNT OR NOT C IF (LTAB(M,12).NE.0) THEN DO 240 IL=1,NGF IF (LTAB(M,1).EQ.IG(IL)) THEN W(1)=WX(IL) W(2)=WY(IL) W(3)=WZ(IL) C FORM (MISCLOSURE VECTOR OF THE CONTROL POINT OBSERVATIONS) IF (LTAB(M,12).EQ.1) THEN WC(1)=XJ-XG(IL) WC(2)=YJ-YG(IL) WC(3)=ZJ-ZG(IL) ELSEIF (LTAB(M,12).EQ.2) THEN WC(1)=XJ-XG(IL) WC(2)=YJ-YG(IL) WC(3)=0.D0 ELSEIF (LTAB(M,12).EQ.3) THEN WC(1)=0.D0 WC(2)=0.D0 WC(3)=ZJ-ZG(IL) ENDIF C COMPLETE AND WITH THE CONTRIBUTION OF THE CONTROL PNT DO 220 JR=1,3 C(JR,JR)=C(JR,JR)+W(JR) AW(JR)=AW(JR)+W(JR)*WC(JR) 220 CONTINUE GO TO 260 ENDIF 240 CONTINUE ENDIF 260 CONTINUE C = -INV CALL SPIN (C,3,3,DET,IDEXP) C C ALL NECESSARY MATRICES FOR PNT M HAVE BEEN COMPUTED C COMPUTE THE ELEMENTS OF THE SUMMATION ++ C SUM1=P4E0(1)+Q1I0(1)+AW(1) SUM2=P4E0(2)+Q1I0(2)+AW(2) SUM3=P4E0(3)+Q1I0(3)+AW(3) C C CRX,Y,Z = * ( + + ). CORRECTIONS C CRX=C(1,1)*SUM1+C(1,2)*SUM2+C(1,3)*SUM3 CRY=C(2,1)*SUM1+C(2,2)*SUM2+C(2,3)*SUM3 CRZ=C(3,1)*SUM1+C(3,2)*SUM2+C(3,3)*SUM3 C C ADD THE CORRECTIONS TO THE APPRX. OBJECT COORDS C XJ=XJ-CRX YJ=YJ-CRY ZJ=ZJ-CRZ C FORM QVB IF(ITER.EQ.ITRF) THEN IF(PAUX.GT.0.D0)THEN QVB=QVB+PAUX*(CRX**2+CRY**2+CRZ**2) WRITE(IW,*) ' QVB = ',QVB ENDIF ENDIF C STORE THE UPDATED OBJECT COORDS C WRITE(UNIT=IEX2,REC=M,ERR=5000) IDO,XJ,YJ,ZJ WRITE(IW,2000) IDO,XJ,YJ,ZJ,CRX,CRY,CRZ C 1000 CONTINUE 2000 FORMAT(' ',/,3X,I10,5X,F12.3,3X,F12.3,3X,F12.3,5X,F8.3,3X, @ F8.3,3X,F8.3) RETURN 5000 WRITE(IW,*) ' **** ERROR RELATED TO FILE ****' RETURN 6000 WRITE(IW,*) ' **** ERROR RELATED TO FILE ****' RETURN END C C ******************************************************************* SUBROUTINE PHORES (W1,QVP,QVI) C C COMPUTES RESIDUALS OF THE PHOTOGRAMMETRIC OBSERVATIONS C * COMPUTES QUADRATIC FORMS OF PHOTO- AND IO RESIDUALS * C * WRITTEN BY: COSTAS ARMENAKIS ON SEPTEMBER 1986 * C ******************************************************************* C IMPLICIT REAL*8 (A-H,O-Z) REAL*8 WI(10) C COMMON /DEV/ IR,IW,IEX1,IEX2,IBYT1,IBYT2 COMMON /ARR/ XL(50),YL(50),DID(100),XG(100),YG(100),ZG(100), @ WX(100),WY(100),WZ(100),XK(100),YK(100),ZK(100), @ IPT(50),IG(100),IK(100), @ INDEX(100),IAPN(450),IIPH(450),LTAB(1000,12) COMMON /CAM/ OE(6),OI(10),OIUD(10),WOI(10),R(3,3) COMMON /VAR/ PX,PY,PAUX,NPHOT,JFP,NGF,ITRF,ITRI,ITER,NC,LDEF C QVI = QUADRATIC FORM OF IO RESIDUALS QVI=0.D0 XSUM=0.D0 YSUM=0.D0 XSUM2=0.D0 YSUM2=0.D0 NUM=0 C FOR EACH PHOTOGRAPH C WRITE(IW,1200) DO 1000 I=1,NPHOT READ (UNIT=IEX1,REC=I,ERR=5000) IPH, @ JPH, @ (OIUD(IL),IL=1,10), @ (OI(IL),IL=1,10), @ (WOI(IL),IL=1,10), @ (OE(IL),IL=1,6), @ IPNF, @ (IPT(IL),XL(IL),YL(IL),IL=1,IPNF) CALL ROTMAT (OE(4),OE(5),OE(6),SOM,SPH,SKP,COM,CPH,CKP,R) C WRITE(IW,1500) IPH C FOR EACH PHOTO-POINT ON PHOTO 'I' DO 500 K=1,IPNF DO 10 L=1,JFP C GET FINAL OBJECT COORDINATES FOR POINT 'K' READ(UNIT=IEX2,REC=L,ERR=6000) IDO,XJ,YJ,ZJ IF (IPT(K).EQ.IDO) THEN DX=XJ-OE(1) DY=YJ-OE(2) DZ=ZJ-OE(3) GO TO 12 ENDIF 10 CONTINUE 12 CONTINUE U1=R(1,1)*DX+R(1,2)*DY+R(1,3)*DZ U2=R(2,1)*DX+R(2,2)*DY+R(2,3)*DZ U3=R(3,1)*DX+R(3,2)*DY+R(3,3)*DZ C SX=XL(K)-OIUD(1) SY=YL(K)-OIUD(2) SXY=SX*SY SX2=SX**2 SY2=SY**2 S2=SX2+SY2 S4=S2**2 S6=S2**3 B=OIUD(6)*S2+OIUD(7)*S4+OIUD(8)*S6 C DX=OIUD(4)*SX+OIUD(5)*SY+B*SX+2.D0*OIUD(9)*SXY+ @ OIUD(10)*(S2+2.D0*SX2) DY=-OIUD(4)*SY+OIUD(5)*SX+B*SY+ @ OIUD(9)*(S2+2.D0*SY2)+2.D0*OIUD(10)*SXY T=OIUD(3)/U3 VXP=(OIUD(1)+DX-T*U1-XL(K))*1000.D0 VYP=(OIUD(2)+DY-T*U2-YL(K))*1000.D0 XSUM=XSUM+VXP YSUM=YSUM+VYP XSUM2=XSUM2+VXP**2 YSUM2=YSUM2+VYP**2 NUM=NUM+1 WRITE(IW,2000) IPT(K),VXP,VYP 500 CONTINUE C COMPUTE RESIDUALS OF IO OBS IF (ITRI.LE.ITRF) THEN WRITE(IW,1150) C CALL WIMAT (OIUD,OI,WI) C COMPUTE QUADRATIC FORM OF IO RESIDUALS DO 600 IL=1,10 QVI=QVI+WI(IL)**2*WOI(IL) WRITE(IW,1160) WI(IL) 600 CONTINUE ENDIF 1000 CONTINUE C COMPUTE THE MEAN AND STANDARD DEVIATION OF THE RESIDUALS TOT=DFLOAT(NUM) XAVG=XSUM/TOT YAVG=YSUM/TOT DEF=TOT-1.D0 XVAR=(XSUM2-XAVG*XSUM)/DEF YVAR=(YSUM2-YAVG*YSUM)/DEF XSTD=DSQRT(XVAR) YSTD=DSQRT(YVAR) WRITE(IW,3000) XAVG,XSTD WRITE(IW,3500) YAVG,YSTD C COMPUTE QUADRATIC FORM OF PHOTO-RESIDUALS QVP=(PX*XSUM2+PY*YSUM2)/(1000.D0**2) WRITE(IW,*) ' QVP =',QVP 1150 FORMAT(' ',/,50X,'RESIDUALS OF IO OBSERVATIONS',/) 1160 FORMAT(' ',10X,F20.10) 1200 FORMAT(' ',/,/,/,50X,'PHOTO-RESIDUALS (IN MM)',/,/,/) 1500 FORMAT(' ',/,/,40X,'PHOTOGRAPH # ',I10,/,' ',15X,'PNT',10X,'VXP', @ 12X,'VYP') 2000 FORMAT(' ',15X,I10,2F10.3) 3000 FORMAT(' ',/,/,20X,'MEAN OF VXP ==> ',F10.3,5X,'STD OF VXP ==>', @ ' +/-',F10.3,/) 3500 FORMAT(' ',20X,'MEAN OF VYP ==> ',F10.3,5X,'STD OF VYP ==> +/-', @ F10.3,/) RETURN 5000 WRITE(IW,5001) 5001 FORMAT('0',/,/,15X,'ERROR RELATED TO FILE ') RETURN 6000 WRITE(IW,6001) 6001 FORMAT('0',/,/,15X,'ERROR RELATED TO FILE ') RETURN END C C ******************************************************************* SUBROUTINE CEKRES C C COMPUTES RESIDUALS AT THE CHECK POINTS C * WRITTEN BY: COSTAS ARMENAKIS ON SEPTEMBER 1986 * C ******************************************************************* C IMPLICIT REAL*8 (A-H,O-Z) C COMMON /DEV/ IR,IW,IEX1,IEX2,IBYT1,IBYT2 COMMON /ARR/ XL(50),YL(50),DID(100),XG(100),YG(100),ZG(100), @ WX(100),WY(100),WZ(100),XK(100),YK(100),ZK(100), @ IPT(50),IG(100),IK(100), @ INDEX(100),IAPN(450),IIPH(450),LTAB(1000,12) COMMON /VAR/ PX,PY,PAUX,NPHOT,JFP,NGF,ITRF,ITRI,ITER,NC,LDEF C XSUM=0.D0 YSUM=0.D0 ZSUM=0.D0 XSUM2=0.D0 YSUM2=0.D0 ZSUM2=0.D0 NUM=0 C FOR EACH CHECK POINT C WRITE(IW,500) DO 100 II=1,NC DO 10 JJ=1,JFP C GET FINAL OBJECT COORDINATES FOR POINT READ(UNIT=IEX2,REC=JJ,ERR=6000) IDO,XJ,YJ,ZJ IF (IK(II).EQ.IDO) THEN RXC=XJ-XK(II) RYC=YJ-YK(II) RZC=ZJ-ZK(II) XSUM=XSUM+RXC YSUM=YSUM+RYC ZSUM=ZSUM+RZC XSUM2=XSUM2+RXC**2 YSUM2=YSUM2+RYC**2 ZSUM2=ZSUM2+RZC**2 NUM=NUM+1 WRITE(IW,550) IK(II),RXC,RYC,RZC GO TO 12 ENDIF 10 CONTINUE 12 CONTINUE 100 CONTINUE C COMPUTE THE MEAN AND STANDARD DEVIATION OF THE RESIDUALS TOT=DFLOAT(NUM) XAVG=XSUM/TOT YAVG=YSUM/TOT ZAVG=ZSUM/TOT DEF=TOT-1.D0 XVAR=(XSUM2-XAVG*XSUM)/DEF YVAR=(YSUM2-YAVG*YSUM)/DEF ZVAR=(ZSUM2-ZAVG*ZSUM)/DEF XSTD=DSQRT(XVAR) YSTD=DSQRT(YVAR) ZSTD=DSQRT(ZVAR) C WRITE(IW,3000) XAVG,XSTD WRITE(IW,3500) YAVG,YSTD WRITE(IW,4000) ZAVG,ZSTD 500 FORMAT(' ',/,/,/,40X,'CHECK-POINT RESIDUALS (IN METRES)',/,/,/, @ ' ',20X,'PNT',10X,'RXC',10X,'RYC',10X,'RZC',/,/) 550 FORMAT(' ',15X,I10,3F10.3) 3000 FORMAT(' ',/,/,20X,'MEAN OF RXC ==> ',F10.3,5X,'STD OF RXC ==>', @ ' +/-',F10.3,/) 3500 FORMAT(' ',20X,'MEAN OF RYC ==> ',F10.3,5X,'STD OF RYC ==> +/-', @ F10.3,/) 4000 FORMAT(' ',20X,'MEAN OF RZC ==> ',F10.3,5X,'STD OF RZC ==> +/-', @ F10.3,/) RETURN 6000 WRITE(IW,6001) 6001 FORMAT('0',/,/,15X,'ERROR RELATED TO FILE ') RETURN END C C ********************************************************************** SUBROUTINE COVAR (JU,NDIM,N,SOC,A,AT,C,C1,E,P4,P,Q1,E2T,ET,E2,P1, @ Q2,PT,P2) C C * COMPUTES THE VARIANCE-COVARIANCE MATRIX OF THE OBJECT POINTS * C * WRITTEN BY: COSTAS ARMENAKIS ON SEPTEMBER 1986 * C ********************************************************************** C IMPLICIT REAL*8 (A-H,O-Z) REAL*8 N(NDIM),W(3),WC(3),A(2,3),AT(3,2),C(3,3),C1(3,3),E(2,6), @ P4(3,6),P(2,10),Q1(3,10),E1(6,6),E2T(3,6),B(3,6),ET(6,2), @ E2(6,3),H1(3,3),H(3,3),S(3,3),CT(3,3),CS(3,3),CM(3,3), @ P1(10,6),F(3,6),D(6,10),SS(3,10),QD(3,10),Q2(10,10), @ Z(3,10),PT(10,2),P2(10,3) C COMMON /ARR/ XL(50),YL(50),DID(100),XG(100),YG(100),ZG(100), @ WX(100),WY(100),WZ(100),XK(100),YK(100),ZK(100), @ IPT(50),IG(100),IK(100), @ INDEX(100),IAPN(450),IIPH(450),LTAB(1000,12) COMMON /CAM/ OE(6),OI(10),OIUD(10),WOI(10),R(3,3) COMMON /DEV/ IR,IW,IEX1,IEX2,IBYT1,IBYT2 COMMON /VAR/ PX,PY,PAUX,NPHOT,JFP,NGF,ITRF,ITRI,ITER,NC,LDEF C FORM INVERSE OF MATRIX C FORM N <== L-INV DO 101 IL=1,JU DO 201 JR=1,IL C THE DIAGONAL ELEMENTS OF ARE ALREADY IN RECIPROCAL FORM C FROM ROUTINE IF (JR.LT.IL) THEN SUM=0.D0 KH=IL-1 DO 301 K1=JR,KH LOC1=(IL*(IL-1)/2)+K1 LOC2=(K1*(K1-1)/2)+JR SUM=SUM-N(LOC1)*N(LOC2) 301 CONTINUE LOC3=(IL*(IL-1)/2)+JR LOC4=(IL*(IL-1)/2)+IL N(LOC3)=SUM*N(LOC4) ENDIF 201 CONTINUE 101 CONTINUE C FORM N <== (L-INV)-TRAN * L-INV DO 700 IL=1,JU DO 600 JR=1,IL SUM=0.D0 DO 500 K1=IL,JU LOC1=(K1*(K1-1)/2)+IL LOC2=(K1*(K1-1)/2)+JR SUM=SUM+N(LOC1)*N(LOC2) 500 CONTINUE LOC3=(IL*(IL-1)/2)+JR N(LOC3)=SUM 600 CONTINUE 700 CONTINUE C FOR EACH POINT K IN DO 8000 K=1,JFP C INITIALIZE MATRICES TO ZERO DO 20 IL=1,3 DO 10 JR=1,3 H(IL,JR)=0.D0 C(IL,JR)=0.D0 10 CONTINUE 20 CONTINUE C GET THE APPRX. COORDS OF PNT K READ(UNIT=IEX2,REC=K,ERR=8500) IDO,XJ,YJ,ZJ C GET PHOTOS ON WHICH PNT K APPEARS DO 7000 L=1,LTAB(K,11) DO 35 IL=1,3 DO 30 JR=1,6 B(IL,JR)=0.D0 F(IL,JR)=0.D0 30 CONTINUE 35 CONTINUE IF (ITRI.LE.ITRF) THEN DO 37 IL=1,3 DO 36 JR=1,10 SS(IL,JR)=0.D0 Z(IL,JR)=0.D0 36 CONTINUE 37 CONTINUE ENDIF DO 6000 J=1,LTAB(K,11) C READ (UNIT=IEX1,REC=LTAB(K,J+1),ERR=9000) IPH, @ JPH, @ (OIUD(IL),IL=1,10), @ (OI(IL),IL=1,10), @ (WOI(IL),IL=1,10), @ (OE(IL),IL=1,6), @ IPNF, @ (IPT(IL),XL(IL),YL(IL),IL=1,IPNF) CALL ROTMAT (OE(4),OE(5),OE(6),SOM,SPH,SKP,COM,CPH,CKP,R) C DX=XJ-OE(1) DY=YJ-OE(2) DZ=ZJ-OE(3) U1=R(1,1)*DX+R(1,2)*DY+R(1,3)*DZ U2=R(2,1)*DX+R(2,2)*DY+R(2,3)*DZ U3=R(3,1)*DX+R(3,2)*DY+R(3,3)*DZ T=U3**2 C FORM MATRIX (PART.DERIV. WRT EO ELEMENTS) CALL EMAT (OIUD(3),U1,U2,U3,T,DX,DY,DZ,R,SOM,SPH,SKP, @ COM,CPH,CKP,E) C FORM MATRIX (PART.DERIV. WRT OBJECT COORDS) CALL AMAT (OIUD(3),U1,U2,U3,T,DX,DY,DZ,R,A) C = - TRANSP CALL TRNSD (AT,3,A,2,2,3) C = * PX,PY DO 80 IL=1,3 AT(IL,1)=AT(IL,1)*PX AT(IL,2)=AT(IL,2)*PY 80 CONTINUE C = * CALL MMULD (P4,3,AT,3,E,2,3,2,6) C EXTRACT Q11(J,L) FROM INVERSE & STORE IT IN JEND=J*6 JBEG=JEND-5 LEND=L*6 LBEG=LEND-5 JJ=0 LL=0 IF (J.GE.L) THEN DO 100 IL=JBEG,JEND JJ=JJ+1 DO 85 JR=LBEG,LEND LL=LL+1 IF (IL.GE.JR) THEN LOC=(IL*(IL-1)/2)+JR E1(JJ,LL)=N(LOC) IF (J.EQ.L) THEN E1(LL,JJ)=E1(JJ,LL) ENDIF ENDIF 85 CONTINUE LL=0 100 CONTINUE ELSE DO 120 IL=LBEG,LEND JJ=JJ+1 DO 110 JR=JBEG,JEND LL=LL+1 IF (IL.GE.JR) THEN LOC=(IL*(IL-1)/2)+JR E1(LL,JJ)=N(LOC) ENDIF 110 CONTINUE LL=0 120 CONTINUE ENDIF C = * CALL MMULD (E2T,3,P4,3,E1,6,3,6,6) C FORM B SUBMATRICES DO 140 IL=1,3 DO 130 JR=1,6 B(IL,JR)=B(IL,JR)+E2T(IL,JR) 130 CONTINUE 140 CONTINUE C CONSIDER IO ELEMENTS IF (ITRI.LE.ITRF) THEN C EXTRACT PHOTO-COORDS OF PNT LTAB(K,1) ON PHOTO LTAB(K,J+1) DO 50 IL=1,IPNF IF (LTAB(K,1).EQ.IPT(IL)) THEN XPHOT=XL(IL) YPHOT=YL(IL) GO TO 60 ENDIF 50 CONTINUE 60 CONTINUE C FORM

MATRIX (PART.DERIV. WRT IO ELEMENTS) CALL PMAT (OIUD,U1,U2,U3,XPHOT,YPHOT,P) C = *

CALL MMULD (Q1,3,AT,3,P,2,3,2,10) C EXTRACT FROM INVERSE & STORE IT IN J6=6*NPHOT JBEG=J6+1+10*(J-1) JEND=J6+10*J LEND=6*L LBEG=LEND-5 JJ=0 LL=0 DO 410 IL=JBEG,JEND JJ=JJ+1 DO 400 JR=LBEG,LEND LL=LL+1 LOC=(IL*(IL-1)/2)+JR P1(JJ,LL)=N(LOC) 400 CONTINUE LL=0 410 CONTINUE C = * CALL MMULD (E2T,3,Q1,3,P1,10,3,10,6) C FORM MATRICES DO 430 IL=1,3 DO 420 JR=1,6 F(IL,JR)=F(IL,JR)+E2T(IL,JR) 420 CONTINUE 430 CONTINUE C EXTRACT -TRANSP FROM INVERSE & STORE IT IN JJ=0 LL=0 LBEG=J6+1+10*(L-1) LEND=J6+10*L JEND=6*J JBEG=JEND-5 DO 450 IL=LBEG,LEND JJ=JJ+1 DO 440 JR=JBEG,JEND LL=LL+1 LOC=(IL*(IL-1)/2)+JR D(LL,JJ)=N(LOC) 440 CONTINUE LL=0 450 CONTINUE C = * CALL MMULD (QD,3,P4,3,D,6,3,6,10) C FORM SUBMATRICES DO 470 IL=1,3 DO 460 JR=1,10 SS(IL,JR)=SS(IL,JR)+QD(IL,JR) 460 CONTINUE 470 CONTINUE C EXTRACT FROM INVERSE & STORE IT IN JBEG=J6+1+10*(J-1) JEND=J6+10*J LBEG=J6+1+10*(L-1) LEND=J6+10*L JJ=0 LL=0 IF (J.GE.L) THEN DO 490 IL=JBEG,JEND JJ=JJ+1 DO 480 JR=LBEG,LEND LL=LL+1 IF (IL.GE.JR) THEN LOC=(IL*(IL-1)/2)+JR Q2(JJ,LL)=N(LOC) IF (J.EQ.L) THEN Q2(LL,JJ)=Q2(JJ,LL) ENDIF ENDIF 480 CONTINUE LL=0 490 CONTINUE ELSE DO 510 IL=LBEG,LEND JJ=JJ+1 DO 505 JR=JBEG,JEND LL=LL+1 IF (IL.GE.JR) THEN LOC=(IL*(IL-1)/2)+JR Q2(LL,JJ)=N(LOC) ENDIF 505 CONTINUE LL=0 510 CONTINUE ENDIF C = * CALL MMULD (QD,3,Q1,3,Q2,10,3,10,10) C FORM SUBMATRICES DO 530 IL=1,3 DO 520 JR=1,10 Z(IL,JR)=Z(IL,JR)+QD(IL,JR) 520 CONTINUE 530 CONTINUE ENDIF C 6000 CONTINUE C READ (UNIT=IEX1,REC=LTAB(K,L+1),ERR=9000) IPH, @ JPH, @ (OIUD(IL),IL=1,10), @ (OI(IL),IL=1,10), @ (WOI(IL),IL=1,10), @ (OE(IL),IL=1,6), @ IPNF, @ (IPT(IL),XL(IL),YL(IL),IL=1,IPNF) CALL ROTMAT (OE(4),OE(5),OE(6),SOM,SPH,SKP,COM,CPH,CKP,R) C DX=XJ-OE(1) DY=YJ-OE(2) DZ=ZJ-OE(3) U1=R(1,1)*DX+R(1,2)*DY+R(1,3)*DZ U2=R(2,1)*DX+R(2,2)*DY+R(2,3)*DZ U3=R(3,1)*DX+R(3,2)*DY+R(3,3)*DZ T=U3**2 C FORM MATRIX (PART.DERIV. WRT EO ELEMENTS) CALL EMAT (OIUD(3),U1,U2,U3,T,DX,DY,DZ,R,SOM,SPH,SKP, @ COM,CPH,CKP,E) C FORM MATRIX (PART.DERIV. WRT OBJECT COORDS) CALL AMAT (OIUD(3),U1,U2,U3,T,DX,DY,DZ,R,A) C = - TRANSP CALL TRNSD (ET,6,E,2,2,6) C = * PX,PY DO 150 IL=1,3 ET(IL,1)=ET(IL,1)*PX ET(IL,2)=ET(IL,2)*PY 150 CONTINUE C = * CALL MMULD (E2,6,ET,6,A,2,6,2,3) C

= * CALL MMULD (H1,3,B,3,E2,6,3,6,3) C UPDATE USING THE 1ST TERM DO 170 IL=1,3 DO 160 JR=1,3 H(IL,JR)=H(IL,JR)+H1(IL,JR) 160 CONTINUE 170 CONTINUE IF (ITRI.LE.ITRF) THEN C

= * CALL MMULD (H1,3,F,3,E2,6,3,6,3) C UPDATE USING 2ND TERM DO 550 IL=1,3 DO 540 JR=1,3 H(IL,JR)=H(IL,JR)+H1(IL,JR) 540 CONTINUE 550 CONTINUE C EXTRACT PHOTO-COORDS OF PNT LTAB(K,1) ON PHOTO LTAB(K,L+1) DO 560 IL=1,IPNF IF (LTAB(K,1).EQ.IPT(IL)) THEN XPHOT=XL(IL) YPHOT=YL(IL) GO TO 570 ENDIF 560 CONTINUE 570 CONTINUE C FORM

MATRIX (PART.DERIV. WRT IO ELEMENTS) CALL PMAT (OIUD,U1,U2,U3,XPHOT,YPHOT,P) C =

-TRANSP CALL TRNSD (PT,10,P,2,2,10) C = * PX,PY DO 580 IL=1,10 PT(IL,1)=PT(IL,1)*PX PT(IL,2)=PT(IL,2)*PY 580 CONTINUE C = * CALL MMULD (P2,10,PT,10,A,2,10,2,3) C

= * CALL MMULD (H1,3,SS,3,P2,10,3,10,3) C UPDATE USING 3RD TERM DO 605 IL=1,3 DO 590 JR=1,3 H(IL,JR)=H(IL,JR)+H1(IL,JR) 590 CONTINUE 605 CONTINUE C

= * CALL MMULD (H1,3,Z,3,P2,10,3,10,3) C UPDATE USING 4TH TERM DO 620 IL=1,3 DO 610 JR=1,3 H(IL,JR)=H(IL,JR)+H1(IL,JR) 610 CONTINUE 620 CONTINUE ENDIF C = - TRANSP CALL TRNSD (AT,3,A,2,2,3) C = * PX,PY DO 180 IL=1,3 AT(IL,1)=AT(IL,1)*PX AT(IL,2)=AT(IL,2)*PY 180 CONTINUE C = * CALL MMULD (C1,3,AT,3,A,2,3,2,3) C UPDATE DO 210 IL=1,3 DO 190 JR=1,3 C(IL,JR)=C(IL,JR)+C1(IL,JR) 190 CONTINUE 210 CONTINUE C 7000 CONTINUE C CONSIDER THE WEIGHT FACTOR OF THE APPRX. OBJECT COORDS IF (PAUX.NE.0.D0) THEN DO 200 IL=1,3 C(IL,IL)=C(IL,IL)+PAUX 200 CONTINUE ENDIF C CHECK IF PNT K=LTAB(K,1) IS A CONTROL PNT OR NOT IF (LTAB(K,12).NE.0) THEN DO 240 IL=1,NGF IF (LTAB(K,1).EQ.IG(IL)) THEN W(1)=WX(IL) W(2)=WY(IL) W(3)=WZ(IL) C COMPLETE WITH THE CONTRIBUTION OF THE CONTROL PNT DO 220 JR=1,3 C(JR,JR)=C(JR,JR)+W(JR) 220 CONTINUE GO TO 260 ENDIF 240 CONTINUE ENDIF 260 CONTINUE C FORM = (H+ATPA+PC) = (H+C) DO 280 IL=1,3 DO 270 JR=1,3 S(IL,JR)=H(IL,JR)+C(IL,JR) 270 CONTINUE 280 CONTINUE C = -INVER CALL SPIN (C,3,3,DET,IDEXP) C = * CALL MMULD (CS,3,C,3,S,3,3,3,3) C = -TRANSP CALL TRNSD (CT,3,C,3,3,3) C = * ===> RELATIVE VARIANCE-COVARIANCE MATRIX CALL MMULD (CM,3,CS,3,CT,3,3,3,3) C ALL NECESSARY MATRICES FOR PNT K HAVE BEEN COMPUTED DO 300 IL=1,3 DO 290 JR=1,3 CM(IL,JR)=SOC*CM(IL,JR) 290 CONTINUE 300 CONTINUE WRITE(IW,*) 'VARIANCE-COVARIANCE MATRIX OF PNT #',LTAB(K,1) CALL MOUTD (CM,3,3,3) XX = DSQRT (CM(1,1)) YY = DSQRT (CM(2,2)) ZZ = DSQRT (CM(3,3)) WRITE(IW,8010) XX,YY,ZZ C STORE PNT ID, X- Y- Z- COORDS,ELEMENTS OF VAR-COV MATRIX(SYMM. MODE) WRITE(UNIT=IEX2,REC=K,ERR=8500) @ IDO,XJ,YJ,ZJ,((CM(IL,JR),JR=1,IL),IL=1,3) 8000 CONTINUE 8010 FORMAT(' ',/,10X,'SX = ',F8.4,2X,'SY = ',F8.4,2X,'SZ = ',F8.4) RETURN 8500 WRITE(IW,*) ' **** ERROR RELATED TO FILE ****' RETURN 9000 WRITE(IW,*) ' **** ERROR RELATED TO FILE ****' RETURN END C C ********************************************************************** SUBROUTINE CEKCON (W,WC,QVP,QVI,QVB,SOC) C C * COMPUTES RESIDUALS AT CONTROL PNTS AND THEIR QUADRATIC FORM * C * COMPUTES A-POSTERIORI VARIANCE FACTOR * C * WRITTEN BY: COSTAS ARMENAKIS ON SEPTEMBER 1986 * C ********************************************************************** C IMPLICIT REAL*8 (A-H,O-Z) REAL*8 W(3),WC(3) C COMMON /ARR/ XL(50),YL(50),DID(100),XG(100),YG(100),ZG(100), @ WX(100),WY(100),WZ(100),XK(100),YK(100),ZK(100), @ IPT(50),IG(100),IK(100), @ INDEX(100),IAPN(450),IIPH(450),LTAB(1000,12) COMMON /DEV/ IR,IW,IEX1,IEX2,IBYT1,IBYT2 COMMON /VAR/ PX,PY,PAUX,NPHOT,JFP,NGF,ITRF,ITRI,ITER,NC,LDEF COMMON /HV/ JF,JH,JV C C QVC = QUADRATIC FORM OF CONTROL PNT RESIDUALS WRITE(IW,500) QVC=0.D0 XSUM=0.D0 YSUM=0.D0 ZSUM=0.D0 XSUM2=0.D0 YSUM2=0.D0 ZSUM2=0.D0 C FOR EACH POINT K IN DO 8000 K=1,JFP C GET THE FINAL COORDS OF PNT K READ(UNIT=IEX2,REC=K,ERR=8500) IDO,XJ,YJ,ZJ C CHECK IF PNT K=LTAB(K,1) IS A CONTROL PNT OR NOT IF (LTAB(K,12).NE.0) THEN DO 240 IL=1,NGF IF (LTAB(K,1).EQ.IG(IL)) THEN W(1)=WX(IL) W(2)=WY(IL) W(3)=WZ(IL) C FORM (MISCLOSURE VECTOR OF THE CONTROL POINT OBSERVATIONS) IF (LTAB(K,12).EQ.1) THEN WC(1)=XJ-XG(IL) WC(2)=YJ-YG(IL) WC(3)=ZJ-ZG(IL) XSUM=XSUM+WC(1) YSUM=YSUM+WC(2) ZSUM=ZSUM+WC(3) XSUM2=XSUM2+WC(1)**2 YSUM2=YSUM2+WC(2)**2 ZSUM2=ZSUM2+WC(3)**2 WRITE(IW,8110) LTAB(K,1),WC(1),WC(2),WC(3) ELSEIF (LTAB(K,12).EQ.2) THEN WC(1)=XJ-XG(IL) WC(2)=YJ-YG(IL) WC(3)=0.D0 XSUM=XSUM+WC(1) YSUM=YSUM+WC(2) XSUM2=XSUM2+WC(1)**2 YSUM2=YSUM2+WC(2)**2 WRITE(IW,8120) LTAB(K,1),WC(1),WC(2) ELSEIF (LTAB(K,12).EQ.3) THEN WC(1)=0.D0 WC(2)=0.D0 WC(3)=ZJ-ZG(IL) ZSUM=ZSUM+WC(3) ZSUM2=ZSUM2+WC(3)**2 WRITE(IW,8130) LTAB(K,1),WC(3) ENDIF QVC=QVC+W(1)*WC(1)**2+W(2)*WC(2)**2+W(3)*WC(3)**2 GO TO 7999 ENDIF 240 CONTINUE ENDIF 7999 CONTINUE 8000 CONTINUE WRITE(IW,*) ' QVC =',QVC C SOC = A-POSTERIORI VARIANCE FACTOR SOC=(QVP+QVI+QVC+QVB)/DFLOAT(LDEF) C COMPUTE THE MEAN AND STANDARD DEVIATION OF THE RESIDUALS MHZ=JF+JH MVR=JF+JV TOTH=DFLOAT(MHZ) TOTV=DFLOAT(MVR) XAVG=XSUM/TOTH YAVG=YSUM/TOTH ZAVG=ZSUM/TOTV DEFH=TOTH-1.D0 DEFV=TOTV-1.D0 XVAR=(XSUM2-XAVG*XSUM)/DEFH YVAR=(YSUM2-YAVG*YSUM)/DEFH ZVAR=(ZSUM2-ZAVG*ZSUM)/DEFV XSTD=DSQRT(XVAR) YSTD=DSQRT(YVAR) ZSTD=DSQRT(ZVAR) WRITE(IW,3000) XAVG,XSTD WRITE(IW,3500) YAVG,YSTD WRITE(IW,4000) ZAVG,ZSTD WRITE(IW,*) ' SOC =',SOC 500 FORMAT(' ',/,/,/,40X,'CONTROL-POINT RESIDUALS (IN METRES)', @ /,/,/,' ',20X,'PNT',10X,'DX',10X,'DY',10X,'DZ',/,/) 550 FORMAT(' ',15X,I10,3F10.3) 3000 FORMAT(' ',/,/,20X,'MEAN OF DX ==> ',F10.3,5X,'STD OF DX ==>', @ ' +/-',F10.3,/) 3500 FORMAT(' ',20X,'MEAN OF DY ==> ',F10.3,5X,'STD OF DY ==> +/-', @ F10.3,/) 4000 FORMAT(' ',20X,'MEAN OF DZ ==> ',F10.3,5X,'STD OF DZ ==> +/-', @ F10.3,/) 8110 FORMAT(' ',10X,I10,3F10.3) 8120 FORMAT(' ',10X,I10,2F10.3) 8130 FORMAT(' ',10X,I10,20X,F10.3) RETURN 8500 WRITE(IW,*) ' **** ERROR RELATED TO FILE ****' RETURN END C C**************************************************************** SUBROUTINE SPIN(Q,N,MM,DET,IDEXP) C C WRITTEN BY: R.R. STEEVES C***************************************************************** REAL*8 Q(MM,MM),DET,SUM,DSQRT,DABS,RPART,APART DET=0.D0 DO 4 J=1,N DO 4 I=1,J IF(I.EQ.1) GO TO 2 M=I-1 SUM=0.0D0 DO 1 K=1,M 1 SUM=SUM+Q(K,I)*Q(K,J) Q(I,J)=Q(I,J)-SUM 2 IF(I.EQ.J) GO TO 3 Q(I,J)=Q(I,J)/Q(I,I) GO TO 4 3 CONTINUE DET=DET+DLOG10(Q(I,I)) Q(I,I)=DSQRT(Q(I,I)) 4 CONTINUE IDEXP=DET RPART=DET-IDEXP APART=DABS(RPART) IF(APART.LT.1.D-20) DET=1.D0 IF(APART.LT.1.D-20) GO TO 10 DET=10**RPART 10 CONTINUE DO 7 J=1,N DO 7 I=1,J IF(I.LT.J) GO TO 5 Q(J,J)=1.0D0/Q(J,J) GO TO 7 5 SUM=0.0D0 M=J-1 DO 6 K=I,M 6 SUM=SUM-Q(I,K)*Q(K,J) Q(I,J)=SUM/Q(J,J) 7 CONTINUE DO 9 J=1,N DO 9 I=1,J SUM=0.0D0 DO 8 K=J,N 8 SUM=SUM+Q(I,K)*Q(J,K) Q(I,J)=SUM IF(I.EQ.J) GO TO 9 Q(J,I)=SUM 9 CONTINUE RETURN END C C SUBROUTINE SYMSYS (A,M,N,B,IB,IDGT,D1,D2,IER) C PURPOSE - LINEAR EQUATION SOLUTION - POSITIVE DEFINITE C MATRIX - SYMMETRIC STORAGE MODE - SPACE C ECONOMIZER SOLUTION (IMSL) C DIMENSION A(1),B(IB,1) DOUBLE PRECISION A,B,D1,D2 C FIRST EXECUTABLE STATEMENT C INITIALIZE IER IER = 0 C DECOMPOSE A CALL CHOLDP (A,A,N,D1,D2,IER) IF (IER.NE.0) GO TO 9000 C PERFORM ELIMINATION DO 5 I = 1,M CALL ELIM (A,B(1,I),N,B(1,I)) 5 CONTINUE GO TO 9005 9000 CONTINUE CALL MESAG(IER,6HSYMSYS) 9005 RETURN END C C SUBROUTINE CHOLDP (A,UL,N,D1,D2,IER) C PURPOSE - DECOMPOSITION OF A POSITIVE DEFINITE MATRIX - C SYMMETRIC STORAGE MODE C STORAGE MODE. THE DIAGONAL OF L CONTAINS THE C RECIPROCALS OF THE ACTUAL DIAGONAL ELEMENTS (IMSL) C C DIMENSION A(1),UL(1) DOUBLE PRECISION A,UL,D1,D2,ZERO,ONE,FOUR,SIXTN,SIXTH,X,RN DATA ZERO,ONE,FOUR,SIXTN,SIXTH/ * 0.0D0,1.D0,4.D0,16.D0,.0625D0/ C FIRST EXECUTABLE STATEMENT D1=ONE D2=ZERO RN = ONE/(N*SIXTN) IP = 1 IER=0 DO 45 I = 1,N IQ = IP IR = 1 DO 40 J = 1,I X = A(IP) IF (J .EQ. 1) GO TO 10 DO 5 K=IQ,IP1 X = X - UL(K) * UL(IR) IR = IR+1 5 CONTINUE 10 IF (I.NE.J) GO TO 30 D1 = D1*X IF (A(IP) + X*RN .LE. A(IP)) GO TO 50 15 IF (DABS(D1).LE.ONE) GO TO 20 D1 = D1 * SIXTH D2 = D2 + FOUR GO TO 15 20 IF (DABS(D1) .GE. SIXTH) GO TO 25 D1 = D1 * SIXTN D2 = D2 - FOUR GO TO 20 25 UL(IP) = ONE/DSQRT(X) GO TO 35 30 UL(IP) = X * UL(IR) 35 IP1 = IP IP = IP+1 IR = IR+1 40 CONTINUE 45 CONTINUE GO TO 9005 50 IER = 129 9000 CONTINUE CALL MESAG(IER,6HCHOLDP) 9005 RETURN END C C SUBROUTINE ELIM (A,B,N,X) C PURPOSE - ELIMINATION PART OF THE SOLUTION OF AX=B - C POSITIVE DEFINITE MATRIX - SYMMETRIC C STORAGE MODE (IMSL) C DIMENSION A(1),B(1),X(1) DOUBLE PRECISION A,B,X,T,ZERO DATA ZERO/0.0D0/ C FIRST EXECUTABLE STATEMENT C SOLUTION OF LY = B IP=1 IW = 0 DO 15 I=1,N T=B(I) IM1 = I-1 IF (IW .EQ. 0) GO TO 9 IP=IP+IW-1 DO 5 K=IW,IM1 T = T-A(IP)*X(K) IP=IP+1 5 CONTINUE GO TO 10 9 IF (T .NE. ZERO) IW = I IP = IP+IM1 10 X(I)=T*A(IP) IP=IP+1 15 CONTINUE C SOLUTION OF UX = Y N1 = N+1 DO 30 I = 1,N II = N1-I IP=IP-1 IS=IP IQ=II+1 T=X(II) IF (N.LT.IQ) GO TO 25 KK = N DO 20 K=IQ,N T = T - A(IS) * X(KK) KK = KK-1 IS = IS-KK 20 CONTINUE 25 X(II)=T*A(IS) 30 CONTINUE RETURN END C C SUBROUTINE MESAG (IER,NAME) C PURPOSE - PRINT A MESSAGE REFLECTING AN ERROR CONDITION C IER = I+J WHERE C I = 128 IMPLIES TERMINAL ERROR MESSAGE, C I = 64 IMPLIES WARNING WITH FIX MESSAGE, C I = 32 IMPLIES WARNING MESSAGE. C J = ERROR CODE RELEVANT TO CALLING C ROUTINE. (IMSL) INTEGER IER INTEGER NAME(1) C SPECIFICATIONS FOR LOCAL VARIABLES INTEGER I,IEQ,IEQDF,IOUNIT,LEVEL,LEVOLD,NAMEQ(6), * NAMSET(6),NAMUPK(6),NIN,NMTB DATA NAMSET/1HU,1HE,1HR,1HS,1HE,1HT/ DATA NAMEQ/6*1H / DATA LEVEL/4/,IEQDF/0/,IEQ/1H=/ C UNPACK NAME INTO NAMUPK C FIRST EXECUTABLE STATEMENT CALL UNPAC (NAME,6,NAMUPK,NMTB) C GET OUTPUT UNIT NUMBER CALL UNINUM(1,NIN,IOUNIT) C CHECK IER IF (IER.GT.999) GO TO 25 IF (IER.LT.-32) GO TO 55 IF (IER.LE.128) GO TO 5 IF (LEVEL.LT.1) GO TO 30 C PRINT TERMINAL MESSAGE IF (IEQDF.EQ.1) WRITE(IOUNIT,35) IER,NAMEQ,IEQ,NAMUPK IF (IEQDF.EQ.0) WRITE(IOUNIT,35) IER,NAMUPK GO TO 30 5 IF (IER.LE.64) GO TO 10 IF (LEVEL.LT.2) GO TO 30 C PRINT WARNING WITH FIX MESSAGE IF (IEQDF.EQ.1) WRITE(IOUNIT,40) IER,NAMEQ,IEQ,NAMUPK IF (IEQDF.EQ.0) WRITE(IOUNIT,40) IER,NAMUPK GO TO 30 10 IF (IER.LE.32) GO TO 15 C PRINT WARNING MESSAGE IF (LEVEL.LT.3) GO TO 30 IF (IEQDF.EQ.1) WRITE(IOUNIT,45) IER,NAMEQ,IEQ,NAMUPK IF (IEQDF.EQ.0) WRITE(IOUNIT,45) IER,NAMUPK GO TO 30 15 CONTINUE C CHECK FOR UERSET CALL DO 20 I=1,6 IF (NAMUPK(I).NE.NAMSET(I)) GO TO 25 20 CONTINUE LEVOLD = LEVEL LEVEL = IER IER = LEVOLD IF (LEVEL.LT.0) LEVEL = 4 IF (LEVEL.GT.4) LEVEL = 4 GO TO 30 25 CONTINUE IF (LEVEL.LT.4) GO TO 30 C PRINT NON-DEFINED MESSAGE IF (IEQDF.EQ.1) WRITE(IOUNIT,50) IER,NAMEQ,IEQ,NAMUPK IF (IEQDF.EQ.0) WRITE(IOUNIT,50) IER,NAMUPK 30 IEQDF = 0 RETURN 35 FORMAT(19H *** TERMINAL ERROR,10X,7H(IER = ,I3, 1 20H) FROM ROUTINE ,6A1,A1,6A1) 40 FORMAT(27H *** WARNING WITH FIX ERROR,2X,7H(IER = ,I3, 1 20H) FROM ROUTINE ,6A1,A1,6A1) 45 FORMAT(18H *** WARNING ERROR,11X,7H(IER = ,I3, 1 20H) FROM ROUTINE ,6A1,A1,6A1) 50 FORMAT(20H *** UNDEFINED ERROR,9X,7H(IER = ,I5, 1 20H) FROM ROUTINE ,6A1,A1,6A1) C C SAVE P FOR P = R CASE C P IS THE PAGE NAMUPK C R IS THE ROUTINE NAMUPK 55 IEQDF = 1 DO 60 I=1,6 60 NAMEQ(I) = NAMUPK(I) 65 RETURN END C C SUBROUTINE UNPAC (PACKED,NCHARS,UNPAKD,NCHMTB) C PURPOSE - NUCLEUS CALLED BY IMSL ROUTINES THAT HAVE C CHARACTER STRING ARGUMENTS (IMSL) INTEGER NC,NCHARS,NCHMTB C LOGICAL*1 UNPAKD(1),PACKED(1),LBYTE,LBLANK INTEGER*2 IBYTE,IBLANK EQUIVALENCE (LBYTE,IBYTE) DATA LBLANK /1H / DATA IBYTE /1H / DATA IBLANK /1H / C INITIALIZE NCHMTB NCHMTB = 0 C RETURN IF NCHARS IS LE ZERO IF(NCHARS.LE.0) RETURN C SET NC=NUMBER OF CHARS TO BE DECODED NC = MIN0 (129,NCHARS) NWORDS = NC*4 J = 1 DO 110 I = 1,NWORDS,4 UNPAKD(I) = PACKED(J) UNPAKD(I+1) = LBLANK UNPAKD(I+2) = LBLANK UNPAKD(I+3) = LBLANK 110 J = J+1 C CHECK UNPAKD ARRAY AND SET NCHMTB C BASED ON TRAILING BLANKS FOUND DO 200 N = 1,NWORDS,4 NN = NWORDS - N - 2 LBYTE = UNPAKD(NN) IF(IBYTE .NE. IBLANK) GO TO 210 200 CONTINUE NN = 0 210 NCHMTB = (NN + 3) / 4 RETURN END C C SUBROUTINE UNINUM(IOPT,NIN,NOUT) C PURPOSE - TO RETRIEVE CURRENT VALUES AND TO SET NEW C VALUES FOR INPUT AND OUTPUT UNIT C IDENTIFIERS. (IMSL) C C SPECIFICATIONS FOR ARGUMENTS INTEGER IOPT,NIN,NOUT C SPECIFICATIONS FOR LOCAL VARIABLES INTEGER NIND,NOUTD DATA NIND/5/,NOUTD/6/ C FIRST EXECUTABLE STATEMENT IF (IOPT.EQ.3) GO TO 10 IF (IOPT.EQ.2) GO TO 5 IF (IOPT.NE.1) GO TO 9005 NIN = NIND NOUT = NOUTD GO TO 9005 5 NIND = NIN GO TO 9005 10 NOUTD = NOUT 9005 RETURN END //GO.PHOTDAT DD DSN=&&TEMP1, // UNIT=VIO,DISP=NEW, // DCB=(RECFM=FB,LRECL=2200,BLKSIZE=4400,DSORG=DA), // SPACE=(TRK,(10,10),RLSE) //GO.APRXOBJ DD DSN=&&TEMP2, // UNIT=VIO,DISP=NEW, // DCB=(RECFM=FB,LRECL=80,BLKSIZE=6000,DSORG=DA), // SPACE=(TRK,(10,10),RLSE) //*GO.APRXOBJ DD DSN=A.M63224.APRXOBJ, //* UNIT=DASD,DISP=(NEW,CATLG,DELETE), //* DCB=(RECFM=FB,LRECL=80,BLKSIZE=6000,DSORG=DA), //* SPACE=(TRK,(10,10),RLSE) //*GO.APRXOBJ DD DSN=A.M63224.APRXOBJ,DISP=SHR //GO.SYSIN DD * ** ENTER DATA ********* //