C C C P R O G R A M - S P S R C C ---------------------------------------- C SINGLE PHOTO SPACE RESECTION - (FORTRAN) C ---------------------------------------- C C PROGRAM CALCULATES OMEGA, PHI, KAPPA TO 0.0001 DEGREES C C ALSO CAMERA STATION COORDINATES C C AND TILT, SWING, AZIMUTH OF PHOTO C C FROM 3 TO 15 CONTROL POINTS MAY BE USED C C RESIDUALS AT CONTROLS ARE CALCULATED C C ------------------------------------ C C INPUT FORMATS ARE INDICATED BELOW THUS -- * C C DIMENSION W(2),P(2),C(2),XC(2),YC(2),ZC(2),X(20),Y(20),AX(20), 1AY(20),AZ(20),A(40,27),AT(27,40),E(40),QXX(27,27),ATE(27),UNK(27), 2V(40),AM(3,3),TITLE(3),STA(3) C ID=0 PI=4.*ATAN (1.) R=180./PI C FORMAT C READ NO. OF PHOTOS (1ST CARD) I5 * READ 420,NOP 420 FORMAT(I5) KKL=1 5 CONTINUE PRINT 724 724 FORMAT(1H1,//,26X,30HSPACE RESECTION - PROGRAM SPSR //) C C READ TITLE AND PRINT (2ND CARD) 2A6,A4 * READ 710,TITLE 710 FORMAT(2A6,A4) PRINT 711,TITLE 711 FORMAT(1H0,5X,2A6,A4) C C READ FOCAL LENGTH AND PRINT (3RD CARD) F8.3 * READ 410,AF 410 FORMAT(F8.3) PRINT 800,AF 800 FORMAT(1H0,25X,14HFOCAL LENGTH =,F9.3,3H MM,///) PRINT 713 713 FORMAT(1H0,25X,'PHOTO COORDS',9X,'GROUND COORDINATES'/) PRINT 714 714 FORMAT(1H0,5X,'CONTROL POINTS'8X,'X',9X,'Y',10X,'E',11X, 1 'N',9X,' HT'/) C C READ NO. OF CONTROLS (4TH CARD) I5 * READ 420,NL NK=1 DO 10 I=1,1 C READ APPROX OMEGA, PHI, KAPPA (RADIANS) C AND AIR STATION COORDS (5TH CARD) * C 10 READ 430,W(I),P(I),C(I),XC(I),YC(I),ZC(I) 430 FORMAT(F7.2,2F5.2,F8.0,F10.0,F5.0) 30 CONTINUE DO 40 I=1,NL C READ 3 CARDS FOR EACH CONTROL NAME CARD 2A6,A4 * C PHOTO COORDS F13.3,F10.3 * C GD COORDS F12.2,2F10.2 * READ 440,STA,X(I),Y(I),AX(I),AY(I),AZ(I) 440 FORMAT(2A6,A4/F13.3,F10.3/F12.2,2F10.2) 40 PRINT 715,STA,X(I),Y(I),AX(I),AY(I),AZ(I) 715 FORMAT(1H0,5X,2A6,A4,2F10.3,2F12.2,F10.2/) 50 CONTINUE N=NL L=1 NN=1 460 CONTINUE AM(1,1)=COS (P(L))*COS (C(L)) AM(1,2)=COS (W(L))*SIN (C(L))+SIN (W(L))*SIN (P(L))*COS (C(L)) AM(1,3)=SIN (W(L))*SIN (C(L))-COS (W(L))*SIN (P(L))*COS (C(L)) AM(2,1)=-COS (P(L))*SIN (C(L)) AM(2,2)=COS (W(L))*COS (C(L))-SIN (W(L))*SIN (P(L))*SIN (C(L)) AM(2,3)=SIN (W(L))*COS (C(L))+COS (W(L))*SIN (P(L))*SIN (C(L)) AM(3,1)=SIN (P(L)) AM(3,2)=-SIN (W(L))*COS (P(L)) AM(3,3)=COS (W(L))*COS (P(L)) IF(ID-1)520,530,530 520 CONTINUE 60 DO 150 I = NN,N 70 K=I J=K+NL C -------------------------------------------------- C GENERATE THE A MATRIX (THE MATRIX OF COEFFICIENTS) C XX=AX(I)-XC(L) YY=AY(I)-YC(L) ZZ=AZ(I)-ZC(L) Q=XX*AM(3,1)+YY*AM(3,2)+ZZ*AM(3,3) C ---------------------------------------------------- C EVALUATE PARTIALS FOR -X- EQUATION AT ASSUMED VALUES C A(K,1)=(X(I)*(YY*(-AM(3,3))+ZZ*AM(3,2))+AF*(YY*(-AM(1,3))+ZZ*AM(1, 12)))/Q A(K,2)=(X(I)*(XX*COS (P(L))+YY*SIN (W(L))*SIN (P(L))-ZZ*COS (W(L)) 1*SIN (P(L)))+AF*((-XX*SIN (P(L)))*COS (C(L))+YY*SIN (W(L))*COS (P( 2L))*COS (C(L))-ZZ*COS (W(L))*COS (P(L))*COS (C(L))))/Q A(K,3)=AF*(XX*AM(2,1)+YY*AM(2,2)+ZZ*AM(2,3))/Q A(K,4)=(X(I)*AM(3,1)+AF*AM(1,1))/Q*(-1.) A(K,5)=(X(I)*AM(3,2)+AF*AM(1,2))/Q*(-1.) A(K,6)=(X(I)*AM(3,3)+AF*AM(1,3))/Q*(-1.) C ---------------------------------------------------- C EVALUATE PARTIALS FOR -Y- EQUATION AT ASSUMED VALUES C A(J,1)=(Y(I)*(YY*(-AM(3,3))+ZZ*AM(3,2))+AF*(YY*(-AM(2,3))+ZZ*AM(2, 12)))/Q A(J,2)=(Y(I)*(XX*COS (P(L))+YY*SIN (W(L))*SIN (P(L))-ZZ*COS (W(L)) 1*SIN (P(L)))+AF*(XX*SIN (P(L))*SIN (C(L))-YY*SIN (W(L))*COS (P(L)) 2*SIN (C(L))+ZZ*COS (W(L))*COS (P(L))*SIN (C(L))))/Q A(J,3)=AF*((-XX*AM(1,1))-YY*AM(1,2)-ZZ*AM(1,3))/Q A(J,4)=(Y(I)*AM(3,1)+AF*AM(2,1))/Q*(-1.) A(J,5)=(Y(I)*AM(3,2)+AF*AM(2,2))/Q*(-1.) A(J,6)=(Y(I)*AM(3,3)+AF*AM(2,3))/Q*(-1.) C ----------------------------------------------- C GENERATE THE E MATRIX (THE MATRIX OF CONSTANTS) C E(K)=(X(I)*(XX*AM(3,1)+YY*AM(3,2)+ZZ*AM(3,3))+AF*(XX*AM(1,1)+YY*AM 1(1,2)+ZZ*AM(1,3)))/(-Q) E(J)=(Y(I)*(XX*AM(3,1)+YY*AM(3,2)+ZZ*AM(3,3))+AF*(XX*AM(2,1)+YY*AM 1(2,2)+ZZ*AM(2,3)))/(-Q) 150 CONTINUE M=2*NL N=6 C C ---------------------- C LEAST SQUARES SOLUTION C ---------------------- C C TRANSPOSE A TO AT DO 260 I=1,M DO 260 J=1,N 260 AT(J,I)=A(I,J) C -------------------- C COMPUTE QXX = AT * A C DO 270 I = 1,N DO 270 J=1,N QXX(I,J)=0. DO 270 K = 1,M 270 QXX(I,J)=QXX(I,J)+AT(I,K)*A(K,J) C --------------------- C INVERT THE QXX MATRIX C DO 330 K=1,N DO 290 J=1,N IF(J-K) 280,290,280 280 QXX(K,J)=QXX(K,J)/QXX(K,K) 290 CONTINUE QXX(K,K)=1./QXX(K,K) DO 330 I=1,N IF(I-K) 300,330,300 300 DO 320 J=1,N IF(J-K) 310,320,310 310 QXX(I,J)=QXX(I,J)-QXX(I,K)*QXX(K,J) 320 CONTINUE QXX(I,K)=-QXX(I,K)*QXX(K,K) 330 CONTINUE C -------------------- C COMPUTE ATF = AT * E C DO 340 I = 1,N ATE(I)=0. DO 340 K = 1,M 340 ATE(I)=ATE(I)+AT(I,K)*E(K) C ----------------------- C COMPUTE UNK = QXX * ATE C DO 350 I = 1,N UNK(I)=0. DO 350 K = 1,N 350 UNK(I)=UNK(I)+QXX(I,K)*ATE(K) 380 CONTINUE W(1)=UNK(1)+W(1) P(1)=UNK(2)+P(1) C(1)=UNK(3)+C(1) XC(1)=XC(1)+UNK(4) YC(1)=YC(1)+UNK(5) ZC(1)=ZC(1)+UNK(6) NK=NK+1 IF(ABS(UNK(1))-.00003)701,701,50 701 IF(ABS(UNK(2))-.00003)702,702,50 702 IF(ABS(UNK(3))-.00003)703,703,50 703 IF(ABS(UNK(4))-.3)704,704,50 704 IF(ABS(UNK(5))-.3)705,705,50 705 IF(ABS(UNK(6))-.3)706,706,50 706 ID=1 GO TO 460 530 CONTINUE ID=0 NK=NK-1 PRINT 716,NK 716 FORMAT(1H0/,6X,26HNO. OF ITERATIONS REQUIRED,I5//) PRINT 550 550 FORMAT(1H0,5X,22HTHE ORIENTATION MATRIX/) DO 540 I=1,3 540 PRINT 570,(AM(I,J),J=1,3) 570 FORMAT(1H0,19X,5F15.6//) DO 390 I=1,1 WW=W(I)*R PP=P(I)*R CC=C(I)*R 390 CONTINUE PRINT 719 719 FORMAT(1H0,5X,18HORIENTATION ANGLES/) PRINT 717 717 FORMAT(1H0,20X,'OMEGA (DEGREES)',5X,'PHI (DEGREES)',7X, 1 'KAPPA (DEGREES)') PRINT 718,WW,PP,CC 718 FORMAT(1H0,23X,F10.4,9X,F10.4,11X,F10.4//) TILT=ATAN (SQRT (1.-AM(3,3)**2)/AM(3,3))*R ALPHA=ATAN (ABS (AM(3,1))/ABS (AM(3,2)))*R SWING=ATAN (ABS (AM(1,3))/ABS (AM(2,3)))*R IF(AM(1,3))600,600,601 600 IF(AM(2,3))606,606,604 601 IF(AM(2,3))603,603,605 603 SWING=360.-SWING GO TO 606 604 SWING=180.-SWING GO TO 606 605 SWING=180.+SWING 606 CONTINUE IF(AM(3,1))610,610,611 610 IF(AM(3,2))616,616,614 611 IF(AM(3,2))613,613,615 613 ALPHA=360.-ALPHA GO TO 615 614 ALPHA=180.-ALPHA GO TO 616 615 ALPHA=180.+ALPHA 616 CONTINUE PRINT 720 720 FORMAT(1H0,20X,'TILT (DEGREES)',6X,'SWING (DEGREES)',5X, 1 'AZIMUTH (DEGREES)') PRINT 718,TILT,SWING,ALPHA PRINT 721 721 FORMAT(/,1H0,5X,26HCAMERA STATION COORDINATES/) PRINT 722 722 FORMAT(1H0,29X,1HE,15X,1HN,15X,2HHT,/) PRINT 723,XC(1),YC(1),ZC(1) 723 FORMAT(1H0,F34.2,2F16.2/) C C -------------------------------------------------- C COMPUTE PHOTO COORDINATE RESIDUALS V = A * UNK - E C DO 360 I = 1,M V(I)=0. DO 360 K = 1,N 360 V(I)=V(I)+A(I,K)*UNK(K) DO 370 I= 1,M 370 V(I)=V(I)-E(I) PRINT 500 500 FORMAT(//,1H0,5X,30HRESIDUALS IN PHOTO COORDINATES/) PRINT 450,(V(I),I=1,M) 450 FORMAT(1H0,5X,5(2F8.3,' ,')//,6X,5(2F8.3,' ,')//, 1 6X,5(2F8.3,' ,')//,6X,5(2F8.3,' ,')) KKL=KKL+1 IF(KKL-NOP)5,5,461 461 CONTINUE C C WHEN CONTINUE WILL GO TO NEW PAGE TO PRINT HEADING FOR NEXT PHOTO C END