C C C C THIS VERSION WAS DONE ON 01/31/86 C E. CARLSON C C C SUBROUTINE TOGRID(CM,PRC,EPS,A,B,C,D,FI,XL,FN,X,Y) IMPLICIT REAL*8(A-H,O-Z) SINFI=DSIN(FI) COSFI=DCOS(FI) T=SINFI/COSFI TS=T*T DL=(XL-CM)*COSFI ET=EPS*COSFI**2 V=DSQRT(1.+ET) OM=FI+B*SINFI*COSFI*(1.+C*COSFI**2*(1.+D*COSFI**2)) XP=OM*PRC*A X=XP+(0.5D0*PRC/V)*T*DL**2*(1.+(1./12.)*DL**2* & (5.-TS+ET*(9+4.*ET)+2.*DL**2*(1.+TS*(-1.+.00167*TS)))) -FN Y=(PRC/V)*DL*(1.+(1./6.)*DL**2*(1.-TS+ET+.05*DL**2* & (5.+TS*(-18.+TS)+15.*ET*(1.-4.*TS)))) RETURN C ENTRY TOGEO(CM,PRC,EPS,A,F,G,H,X,Y,FN,FI,XL) OM=(X+FN)/(A*PRC) SINFI=DSIN(OM) COSFI=DCOS(OM) FI=OM+F*SINFI*COSFI*(1.+G*COSFI**2*(1.+H*COSFI**2)) SINFI=DSIN(FI) COSFI=DCOS(FI) ET=EPS*COSFI**2 V=DSQRT(1.+ET) T=SINFI/COSFI TS=T*T Q=(V/ PRC)*Y FI=FI+0.5D0*T*Q**2*(-1.-ET+(1./12.)*Q**2*(5.+3.* & (TS*(1.-ET*(2.+3.*ET))+ET*(2.-ET)) & -.5*Q**2*(4.+3.*TS*(2.+TS)))) DLCOS=Q*(1.-(1./6.)*Q*Q*(1.+TS+TS+ET-.05*Q**2*(5.05+4.*TS & *(7+6.*TS)))) XL=CM+DLCOS/COSFI RETURN END C SUBROUTINE DMSTOR(KX,MX,SX,ANGLE) IMPLICIT REAL*8(A-H,O-Z) COMMON ROSEC ANGLE=((DBLE(FLOAT(KX))*60.D0+DBLE(FLOAT(MX)))*60.D0+SX)/ROSEC RETURN C ENTRY RTODMS(ANGLE,ID,M,S) TRIFLE=5.D-06 DEG=(ANGLE*ROSEC+TRIFLE)/3600.D0 ID=DEG DEG=(DEG-DBLE(FLOAT(ID)))*60.D0 388 M=DEG S=(DEG-DBLE(FLOAT(M)))*60.D0-TRIFLE 390 S=DABS(S) RETURN END C C* TRIMX SUBROUTINE TRIMAT(NUNK,NEQ,STER,T,V,C,IN) IMPLICIT REAL*8(A-H,O-Z) INTEGER Z DIMENSION T(NUNK),V(NUNK),C(1) C DEFINITION OF CALL ARGUMENTS C NUNK NUMBER OF UNKNOWNS C NEQ NUMBER OF EQUATIONS ON DEVICE IN C STER STANDARD ERROR OF UNIT WEIGHT C T ARRAY IN WHICH CORRECTIONS ARE RETURNED C V ARRAY IN WHICH PRECISION TERMS ARE RETURNED C C ARRAY IN WHICH THE INVERSE MATRIX IS RETURNED C IN DEVICE CONTAINING OBSERVATION EQUATION C ** C ARRAY MUST BE DIMENSIONED BY THE CALLING PROGRAM C THE SIZE OF C MUST BE NUNK*(NUNK+1)/2 LOCATIONS. C DIMENSION JA(1200),RA(1200) 1 C ** JA=UNKNOWN NR RA=ITS COEF. IF EQ HAS MORE THAN 40 VALUES C THE DIMENSION STATEMENT MUST BE CHANGED. DATA BLANK/6H / C FUNCTION LOC COMPUTES THE I ROW AND J COL LOCATION IN C MATRIX LOC(I,J,K)=((I-1)*(K-I+2))/2+J C JA ARRAY CONTAINS UNK NRS IN OBSERVATION C RA ARRAY CONTAINS COEF FOR JA UNKS C T ARRAY WILL CONTAIN THE CORRECTIONS C V ARRAY WILL CONTAIN THE PRECISION OF UNKNOWNS C C ARRAY TO BE INVERTED THE SIZE IS EQUAL TO (M*M+M)/2 WHERE C M IS EQUAL TO NUMBER OF UNKNOWNS C NUNK=M=NR UNKNOWNS, NEQ=NR OF OBSERVATIONS, STER=STANDARD ERROR C CL = COMPUTED-OBSERVED, P=WEIGHT REWIND IN 52 M=NUNK Z= (M-1)*2 N=NEQ MM= (M*M+M)/2 100 SB=0.0 DO 50 I=1,MM 50 C(I)=0.0 DO 51 I=1,M T(I)=0.0 51 V(I)=0.0 C * * * C= A TRANS WA * * * * C L VECTOR NAZ=N DO 102 L=1,N READ(IN,END=103)NAL,(JA(I),RA(I),I=1,NAL),CL,P TI448 IF(P.EQ.0.) NAZ=NAZ-1 PCL=P*CL K=NAL-I JAI=JA(I) DO 1 J=1,K IF(JA(J).LT.JA(J+1)) GO TO 1 TEM=RA(J) IT=JA(J) RA(J)=RA(J+1) JA(J)=JA(J+1) RA(J+1)=TEM JA(J+1)=IT 1 CONTINUE DO 201 I=1,NAL IF(JA(I).EQ.0) GO TO 201 JAI=JA(I) DO 203 II=I,NAL IF(JA(II).EQ.0) GO TO 203 JAII=JA(II) K=LOC(JAI,JAII,Z) C(K) =C(K) +(RA(I)*RA(II)*P) 203 CONTINUE T(JAI)=T(JAI)-RA(I)*PCL 201 CONTINUE 102 SB=SB+PCL*CL 103 CONTINUE DO 2 L=1,M K=LOC(L,L,Z) IF(C(K).EQ.0.0) GO TO 3 GO TO 2 3 WRITE(6 ,400) L STER=BLANK 2 CONTINUE IF(STER.EQ.BLANK) GO TO 4 C REWIND IN C(1)= DSQRT(C(1)) DO 31 J=2,M 31 C(J)=C(J)/C(1) DO 32 I=2,M S=0.0 I1=I-1 DO 33 L=1,I1 K=LOC(L,I,Z) 33 S=S+C(K)*C(K) K=LOC(I,I,Z) VAL= DSQRT(C(K)-S) C(K)=VAL IF(I.GE.M) GO TO 32 IP=I+1 DO 35 J=IP,M S=0.0 IM=I-1 DO 36 L=1,IM K=LOC(L,I,Z) KK=LOC(L,J,Z) 36 S=S+C(K)*C(KK) K=LOC(I,J,Z) 35 C(K)=(C(K)-S)/VAL 32 CONTINUE C V(1)=T(1)/C(1) DO 37 I=2,M S=0.0 IP=I-1 DO 38 L=1,IP K=LOC(L,I,Z) 38 S=S+C(K)*V(L) K=LOC(I,I,Z) 37 V(I)=(T(I)-S)/C(K) DO 39 I=1,M K=LOC(I,I,Z) C(K)=1./C(K) IF(I.GE.M) GO TO 39 ID=I+1 DO 41 J=ID,M S=0.0 JP=J-1 DO 42 K=I,JP KA=LOC(I,K,Z) KB=LOC(K,J,Z) 42 S=S+C(KA)*C(KB) KC=LOC(I,J,Z) KE=LOC(J,J,Z) 41 C(KC)=-S/C(KE) 39 CONTINUE TA=0.0 DO 43 I=1,M 43 TA=TA+V(I)*V(I) VNM=NAZ-M +1.D-6 STER= DSQRT((SB-TA)/VNM) DO 44 I=1,M S=0.0 TA=0.0 DO 45 J=I,M K=LOC(I,J,Z) VAL=C(K) S=S+V(J)*VAL 45 TA=TA+VAL*VAL T(I)=S 44 V(I)=STER*DSQRT(TA) DO 6 L=1,NUNK DO 6 I=L,NUNK S=0. DO 7 K=I,NUNK LK=LOC(L,K,Z) IK=LOC(I,K,Z) 7 S=S+C(LK)*C(IK) LI=LOC(L,I,Z) 6 C(LI)=S C CORRECTIONS IN T ARRAY, PRE. UNKS IN V ARRAY 4 RETURN 400 FORMAT(24H THERE IS NO EQ. FOR UNK ,I3) END C SUBROUTINE ROSTRA(X1,Y1,X2,Y2,NBASE,NCOMN,SIGMA,ANS,PREC,Q,V, 1 IUNIT) IMPLICIT REAL*8(A-H,O-Z) DIMENSION X1(1),Y1(1),X2(1),Y2(1) DIMENSION IU(3),C(3),ANS(1),PREC(1),Q(1),V(1) P=1. LOBSEQ=3 DO 10 I=1,NBASE IF(X1(I).EQ.0.) GO TO 10 IU(1)=1 C(1)=1. IU(2)=3 C(2)=-Y1(I) IU(3)=4 C(3)=X1(I) CMO=-X2(I) WRITE(IUNIT) LOBSEQ,(IU(J),C(J),J=1,3),CMO,P IU(1)=2 C(2)=X1(I) C(3)=Y1(I) CMO=-Y2(I) WRITE(IUNIT) LOBSEQ,(IU(J),C(J),J=1,3),CMO,P 10 CONTINUE NUNK=4 NOBS=NCOMN*2 IF(NOBS.GT.2) GO TO 11 NOBS=4 DO 12 I=1,4 IU(I)=0 12 C(I)=0. DO 13 I=3,4 IU(1)=I C(1)=1. CMO=0. P=1.D15 13 WRITE(IUNIT) LOBSEQ,(IU(J),C(J),J=1,3),CMO,P 11 CONTINUE C CALL TRIMAT(NUNK,NOBS,SIGMA,ANS,PREC,Q,IUNIT) REWIND IUNIT DO 20 I=1,NBASE IF(X1(I).EQ.0.) GO TO 20 K=I+I-1 READ(IUNIT)LOBSEQ,(IU(J),C(J),J=1,3),CMO,P V(K)=ANS(IU(1))*C(1)+ANS(IU(2))*C(2)+ANS(IU(3))*C(3)+CMO V(K)=-V(K) K=K+1 READ(IUNIT)LOBSEQ,(IU(J),C(J),J=1,3),CMO,P V(K)=ANS(IU(1))*C(1)+ANS(IU(2))*C(2)+ANS(IU(3))*C(3)+CMO V(K)=-V(K) 20 CONTINUE RETURN END C SUBROUTINE REFINE(X,Y,V,NEWFIX,NCOMN,XX,YY,CX,CY) IMPLICIT REAL*8(A-H,O-Z) DIMENSION X(1),Y(1),V(1), RAY(4),IN(4),VX(4),VY(4) N=4 IF(NCOMN.LT.N) N=NCOMN DO 11 I=1,N RAY(I)=1.D12 11 IN(I)=0 C DO 10 I=1,N DO 10 J=1,NEWFIX IF(X(J).EQ.0.) GO TO 10 SSQ=(X(J)-XX)**2 + (Y(J)-YY)**2 IF(SSQ.LT.1.) SSQ=1.D-6 IF(I.EQ.1) GO TO 15 M=I DO 14 K=1,M IF(SSQ.EQ.RAY(K)) GO TO 10 14 CONTINUE 15 CONTINUE IF(SSQ.GE.RAY(I)) GO TO 10 RAY(I)=SSQ IN(I)=J VX(I)=V(J+J-1) VY(I)=V(J+J) 10 CONTINUE SUMP=0. SUMX=0. SUMY=0. DO 12 I=1,N P=1./RAY(I) SUMP=SUMP+P SUMX=SUMX+VX(I)*P SUMY=SUMY+VY(I)*P 12 CONTINUE CX=SUMX/SUMP CY=SUMY/SUMP XX=XX+CX YY=YY+CY RETURN END C SUBROUTINE PAGE(IPAGE,KOUNT,KIND) NR=51 KOUNT=KOUNT+1 IF(KOUNT.GE.NR) KOUNT=1 IF(KOUNT.NE.1) RETURN IPAGE=IPAGE+1 WRITE(6,100)IPAGE 100 FORMAT(1H1,100X,'PAGE',I4,/) GO TO (1,2,3,4),KIND 1 WRITE(6,10) 10 FORMAT(1X,'COORDINATES OF COMMON STATIONS',//,1X,'STATION', & 27X,'OLD LATITUDE',6X,'OLD LONGITUDE',5X,'NEW LATITUDE',6X, & 'NEW LONGITUDE',/) RETURN 2 WRITE(6,20) 20 FORMAT(1X,'RESIDUALS AT COMMON STATIONS', & T77,'RESIDUALS AT',//,1X,'STATION',26X, & 'NEW LATITUDE',5X,'NEW LONGITUDE', 12X,'N',7X,'E') WRITE(6,25) 25 FORMAT(T74,'(METERS)',T82,'(METERS)',/) RETURN 3 WRITE(6,30) 30 FORMAT(1X,'TRANSFORMED COORDINATES AND THEIR CORRECTIONS', & T111,'RESIDUALS AT',//,1X, & 'STATION',27X,'OLD LATITUDE',6X,'OLD LONGITUDE',5X, & 'NEW LATITUDE', 6X, 'NEW LONGITUDE', 8X,'N',9X,'E') WRITE(6,35) 35 FORMAT(T108,'(METERS)',T118,'(METERS)',/) RETURN 4 WRITE(6,40) 40 FORMAT(1X,'CHANGES AT INTERSECTIONS POINTS ',T50, & 'SHIFT = NAD83 DATUM - NAD27 DATUM ',//) WRITE(6,41) 41 FORMAT(T10,'OLD',T24,'OLD',T36,'SHIFT',T48,'SHIFT',T62, & 'NEW',T80,'NEW',T95,'SHIFT',T106,'SHIFT') WRITE(6,42) 42 FORMAT(4X, ' LATITUDE',5X,'LONGITUDE',7X,'LAT.',8X,'LON.', TI686 & T61,'LATITUDE',T77,'LONGITUDE',T95,'METERS',T106,'METERS') WRITE(6,43) 43 FORMAT(T37,'SEC',T49,'SEC',T97,'LAT',T107,'LONG',/) RETURN END C SUBROUTINE AFFINE(X1,Y1,X2,Y2,NBASE,NCOMN,SIGMA,ANS,PREC,Q,V, 1 IUNIT) IMPLICIT REAL*8(A-H,O-Z) DIMENSION X1(1),Y1(1),X2(1),Y2(1) DIMENSION IU(3),C(3),ANS(1),PREC(1),Q(1),V(1) P=1. LOBSEQ=3 DO 10 I=1,NBASE IF(X1(I).EQ.0.) GO TO 10 IU(1)=1 C(1)=1. IU(2)=2 C(2)=Y1(I) IU(3)=3 C(3)=X1(I) CMO=-X2(I) WRITE(IUNIT) LOBSEQ,(IU(J),C(J),J=1,3),CMO,P IU(1)=4 IU(2)=5 IU(3)=6 C(2)=X1(I) C(3)=Y1(I) CMO=-Y2(I) WRITE(IUNIT) LOBSEQ,(IU(J),C(J),J=1,3),CMO,P 10 CONTINUE NUNK=6 NOBS=NCOMN*2 C CALL TRIMAT(NUNK,NOBS,SIGMA,ANS,PREC,Q,IUNIT) REWIND IUNIT DO 20 I=1,NBASE IF(X1(I).EQ.0.) GO TO 20 K=I+I-1 READ(IUNIT)LOBSEQ,(IU(J),C(J),J=1,3),CMO,P V(K)=ANS(IU(1))*C(1)+ANS(IU(2))*C(2)+ANS(IU(3))*C(3)+CMO V(K)=-V(K) K=K+1 READ(IUNIT)LOBSEQ,(IU(J),C(J),J=1,3),CMO,P V(K)=ANS(IU(1))*C(1)+ANS(IU(2))*C(2)+ANS(IU(3))*C(3)+CMO V(K)=-V(K) 20 CONTINUE RETURN END