C PROGRAMME HORLAT C************************************************************************ C PROGRAMME HORLAT * C LATITUDE BY THE HORREBOW-TALCOTT METHOD. * C * C COMPUTES UP TO 150 STAR PAIRS FOR EACH STATION. NUMBER OF PAIRS C CAN BE INCREASED BY CORRESPONDING CHANGE IN DIMENSION STATEMENT. C * C************************************************************************ C * C ORDER OF INPUT DATA CARDS. * C (1) STATION CARD. * C * C COLS. 1-20 STATION NAME. * C COLS. 21-23 DEGREES OF LATITUDE (NEGATIVE SIGN FOR SOUTH). * C COLS. 25-26 MINUTES OF LATITUDE. * C COLS. 28-32 SECONDS OF LATITUDE (2 DECIMALS). * C COLS. 34-36 HOURS OF LONGITUDE (NEGATIVE SIGN FOR WEST) * C COLS. 38-39 MINUTES OF LONGITUDE. * C COLS. 41-45 SECONDS OF LONGITUDE (2 DECIMALS). * C COLS. 47-51 TEMPERATURE (2 DECIMALS). * C COLS. 53-57 PRESSURE (2 DECIMALS). * C * C (2) INSTRUMENT CARD. * C * C COLS. 1-8 INSTRUMENT TYPE. * C COLS. 10-15 INSTRUMENT NUMBER. * C COLS. 17-28 CHRONOMETER TYPE. * C COLS. 30-35 CHRONOMETER NUMBER. * C COLS. 37-40 HORREBOW BUBBLE NUMBER. * C COLS. 42-46 HORREBOW BUBBLE SENSITIVITY. * C COLS. 48-51 HORREBOW BUBBLE NUMBER. * C COLS. 53-57 HORREBOW BUBBLE SENSITIVITY. * C COLS. 59-65 EQUATORIAL VALUE OF DRUM REV. (3 DECIMALS). C * C (3) DATE CARD. * C * C COLS. 2-3 FIRST DAY. * C COLS. 5-6 SECOND DAY. * C COLS. 8-9 MONTH. * C COLS. 11-14 YEAR. * C COLS. 16-31 OBSERVER'S NAME. * C COLS. 33-48 RECORDER'S NAME. * C COLS. 50-55 X-COORDINATE OF MEAN POLE IN SECONDS (3 DECIMALS). * C COLS. 57-62 Y-COORDINATE OF MEAN POLE IN SECONDS (3 DECIMALS). * C * C (4) APPARENT PLACE CARDS. * C * C COLS. 2-6 STAR NUMBER. * C COLS. 8-12 STAR MAGNITUDE. * C COLS. 14-15 HOURS OF RIGHT ASCENSION. * C COLS. 17-18 MINUTES OF RIGHT ASCENSION. * C COLS. 20-25 SECONDS OF RIGHT ASCENSION (3 DECIMALS). * C COLS. 27-29 DEGREES OF DECLINATION. * C COLS. 31-32 MINUTES OF DECLINATION. * C COLS. 34-38 SECONDS OF DECLINATION (2 DECIMALS). * C COLS. 40-45 DIU. ABR. COR. IN RT. ASC. IN SEC. ARC (3 D'MALS) * C COLS. 47-52 DIU. ABR. COR. IN DECLINA. IN SEC. ARC (3 D'MALS) * C * C (5) BLANK CARD TO INDICATE END OF APPARENT PLACE CARDS. * C * C (6) OBSERVATION CARDS. * C COLS. 1-6 STAR MUNBER (SOUTH STAR IS NEGATIVE) * C COLS. 8-10 MICROMETER (EAST STAR IS NEGATIVE). * C COLS. 12-15 DRUM DIVISIONS (1 DECIMAL). * C COLS. 17-21 NORTH END OF HORREBOW BUBBLE (1 DECIMAL). * C COLS. 23-27 SOUTH END OF HORREBOW BUBBLE (1 DECIMAL). * C COLS. 29-33 NORTH END OF HORREBOW BUBBLE (1 DECIMAL). * C COLS. 35-39 SOUTH END OF HORREBOW BUBBLE (1 DECIMAL). * C COLS. 41-44 OFF-MER. COR. IN SEC. OF TIME (2 DECIMALS). * C COL. 80 '0' FOR UPPER TRANSIT, '1' FOR LOWER TRANSIT. * C * C (7) BLANK CARD TO INDICATE END OF OBSERVATIONS FOR ONE NIGHT. * C * C (8) STAR PAIR CARDS. C * C COLS. 2-6 FIRST STAR OF PAIR. * C COLS. 8-12 SECOND STAR OF PAIR. * C * C (9) BLANK CARD TO INDICATE END OF STAR PAIRS FOR ONE NIGHT. * C * C * C************************************************************************ IMPLICIT REAL * 8(A-H,O-Z) REAL * 8 RA(150),ULAT(150),ALAT(150),V(150),SE(150),DEC(150),SN(15 10),H(150),DM(150),GMIC(150),DMN(150) INTEGER ID(150),MI(150),IREJ(150),STAR(150),REJ1,REJ2,REJ3,LSTAR(150 150),IST1(150), IST2(150),F1,TH1,NST1(150),NST2(150) DIMENSION AMAT(150,2),ATPA(2,2),ATPW(2),XMAT(2),AX(150),VM(150), 1DMM(150),P(150) INTEGER BLANK/' '/ PI=3.14159265358979D0 RHO=206264.80625 TWOPI=PI*2. C C READING STATION CARD-STATION NAME,LATITUDE(IN DEGREES),LONGITUDE C (IN TIME),TEMPERATURE,PRESSURE. C 19 READ(5,11) F1,F2,F3,F4,F5,IPS,MPS,SPS,ILS,MLS,SLS,TEM,PRE C C READING INSTRUMENT CARD-THEODOLITE TYPE AND NUMBER,CLOCK TYPE AND C NUMBER,HORREBOW BUBBLES NUMBER AND SENSITIVITIES,LEV. SENSITIVITY C 91 READ(5,9) TH1,TH2,NUMT,CL1,CL2,CL3,NUMC,IHR1,BUB1,IHR2,BUB2,R 9 FORMAT(2A4,I7,1X,3A4,I7,I5,F6.3,I5,F6.3,F8.3) WRITE(6,53) 53 FORMAT('1',39X,'LATITUDE BY THE HORREBOW-TALCOTT METHOD',///) C=0. DR=0. M=0 N=0 11 FORMAT(5A4,2I3,F6.2,I4,I3,3F6.2) STALAT=(IPS*3600.D0+MPS*60.D0+SPS)/RHO C C READING DATE CARD-DATE,OBSERVER,RECORDER,POLAR MOTION COORDINATES. C 93 READ(5,8) IDAY1,IDAY2,MO,IYEAR,O1,O2,O3,O4,R1,R2,R3,R4,XP,YP 8 FORMAT(3I3,I5,2(1X,4A4),2F7.3) IF(IYEAR.EQ.0) GO TO 101 C C READING APPARENT PLACE CARDS-STAR NO,RA,DECLINATION,DIURNAL C ABERATION CORRECTIONS,MAGNITUDE,CLOCK CORRECTION. C C NCAT=NO OF UPDATED CATALOGUE STARS. NCAT=0 DO 7 I=1,150 READ(5,4) STAR(I),FMAG,IP,MP,SP,IL,ML,SL,DRA,DDEC 4 FORMAT(' ',I5,F6.2,2I3,F7.3,I4,I3,F6.2,2F7.3) IF(STAR(I).EQ.0) GO TO 81 NCAT=NCAT+1 RA(I)= (IP*3600.D0+MP*60.D0+SP)*15/RHO IF(IL.LT.0) GO TO 30 DEC(I)= (IL*3600.D0+ML*60.D0+SL) /RHO GO TO 7 30 DEC(I)= (IL*3600.D0-ML*60.D0-SL) /RHO 7 CONTINUE 81 IF(M.GT.0) GO TO 50 WRITE(6,13) F1,F2,F3,F4,F5,IPS,MPS,SPS,ILS,MLS,SLS,TH1,TH2,NUMT,O1 1,O2,O3,O4,CL1,CL2,CL3,NUMC,R1,R2,R3,R4,R 13 FORMAT(' ',24X,'LOCATION :',5A4,14X,'INSTRUMENTS',/,25X,'LATITUDE 1:',I4,I3,F7.3,/,25X,'LONG. :',I4,I3,F7.3,12X,2A4,' THEODOLITE N 2O ',I8,/,25X,'OBSERVER :',4A4,10X,3A4,' CHRON. NO',I8,/,25X,'RECOR 3DER :',4A4,10X,'DRUM VALUE',F8.3,'SECS.',///) R=R/2.D0 50 WRITE(6,51) IDAY1,IDAY2,MO,IYEAR 51 FORMAT(' ',////////,50X,'DATE ',I6,'-',I2,'/',I2,'/',I4,//) WRITE(6,44) 44 FORMAT(' ',38X,'STAR POS.OCC. MICROMETER LEVEL 1',/,40X,'NO TURNS DIV NORTH SOUTH',//) M=M+1 C C READ STAR OBSERVATION CARDS-STAR NUMBER,MICR.TURNS,MICR,DIVISIONS, C LEVEL NORTH,LEVEL SOUTH,LEVEL NORTH,LEVEL SOUTH,HOUR ANGLE(IN SECS) C INDEX(IF STAR IS OBSERVED AT LOWER TRANSIT,PUNCH '1' IN COL. ,0) C C ISTAR,LSTAR,ARE OBSERVED STARS C LSTAR IS POSITIVE FOR NORTH STARS C MICR IS POSITIVE FOR STARS OBSERVED IN OCULAR WEST. C C KOUNT=NO OF STARS OBSERVED C KOUNT=0 DO 15 I=1,100 READ(5,14) LSTAR(I),MICR,DIV,ZN1,S1,ZN2,S2,H(I),LOWER 14 FORMAT(I6,I4,F5.1,4F6.1,F5.2,I36) KOUNT=KOUNT+1 H(I)=0.5D0*H(I)*15.D0/RHO IF(LSTAR(I).EQ.0) GO TO 70 IF(LSTAR(I).GT.0) GO TO 45 ISTAR =-LSTAR(I) WRITE(6,41) ISTAR 41 FORMAT(' ',37X,I5,3X,'S') GO TO 46 45 ISTAR = LSTAR(I) WRITE(6,40) ISTAR 40 FORMAT(' ',37X,I5,3X,'N') 46 IF (MICR .GT.0) GO TO 47 GMIC(I)=MICR-DIV/100.D0 SN(I)=-(ZN1+S1+ZN2+S2) MIC =-MICR WRITE(6,42) 42 FORMAT('+',49X,'E') GO TO 48 47 WRITE(6,43) 43 FORMAT('+',49X,'W') GMIC(I)=MICR+DIV/100.D0 SN(I)= (ZN1+S1+ZN2+S2) MIC = MICR 48 DO 18 K=1,NCAT IF(STAR(K).EQ.IABS(LSTAR(I))) GO TO 16 18 CONTINUE WRITE(6,17) 17 FORMAT('+',62X,'STAR NOT UPDATED',/) GO TO 15 16 WRITE(6,49) MIC,DIV,ZN1,S1,ZN2,S2 49 FORMAT('+',52X,I3,F9.1,F8.1,F7.1,3X,/,66X,2F7.1,/) 15 CONTINUE C C COMPUTING BUBBLE VALUES C 70 CONTINUE WRITE(6,1000) 1000 FORMAT('1',//////////) WRITE(6,1) 1 FORMAT(' ',19X,'PAIR',2X,'STAR1',3X,'STAR2',6X,'MEAN',7X,'MICROMET 1ER',5X,'LEVEL',2X,'REFRACTION',3X,'UNADJUSTED',/,21X,'NO',19X,'DEC 2LINATION',5X,'(SECS)',7X,'(SECS)',3X,'(SECS)',6X,'LATITUDE',/) IF(BUB1 .EQ.0.) GO TO 98 D=(BUB1 +BUB2 )/8.D0 GO TO 99 98 D=BUB1 /4.D0 99 CONTINUE REJ =180.D0/RHO DO 100 I=1,100 C C READING STAR PAIRS. C READ(5,106) IST1(I),IST2(I) 106 FORMAT(2I6) IF(IST1(I).EQ.0) GO TO 93 102 N=N+1 DO 107 J=1,KOUNT IF(IABS(LSTAR(J)).EQ.IST1(I)) GO TO 60 107 CONTINUE WRITE(6,29) IST1(I) 29 FORMAT(' ',42X,'STAR NO ',I5,1X,'HAS NOT BEEN OBSERVED',//) N=N-1 GO TO 100 60 DO 56 K=1,NCAT IF(STAR(K).EQ.IST1(I)) GO TO 108 56 CONTINUE 103 N=N-1 WRITE(6,109) IST1(I) 109 FORMAT(' ',42X,'STAR NO',I6,2X,'NOT UPDATED',/) GO TO 100 108 DZ1=-DSIN(2.D0*DEC(K))*DSIN(H(J))*DSIN(H(J))*RHO IF(LOWER.NE.1) GO TO 104 DEC(K)=-DEC(K)+PI 104 DEC1 =DEC(K) IF(LOWER.EQ.1) GO TO 110 DZ1=-DZ1 110 GMIC1 =GMIC(J) SN1=SN(J) DO 116 J=1,KOUNT IF(IABS(LSTAR(J)).EQ.IST2(I)) GO TO 61 116 CONTINUE WRITE(6,29) IST2(I) N=N-1 GO TO 100 61 DO 57 K=1,NCAT IF(STAR(K).EQ.IST2(I)) GO TO 112 57 CONTINUE N=N-1 WRITE(6,109) IST2(I) GO TO 100 112 DZ2=-DSIN(2.D0*DEC(K))*DSIN(H(J))*DSIN(H(J))*RHO IF(LOWER.NE.1) GO TO 113 DEC(K)=-DEC(K)+PI 113 DEC2 =DEC(K) IF(LOWER.EQ.1) GO TO 115 DZ2=-DZ2 115 GMIC2=GMIC(J) SN2=SN(J) DM(I)=GMIC1+GMIC2 DL =(SN1+SN2)/2. FLEV=D*DL ZD=DABS(DEC1 - DEC2 )/2. XMIC=DM(I)*R REF =28.95*DSIN(XMIC*2.D0/RHO)/(DCOS(ZD)*DCOS(ZD)) DECMIN =(DEC1 +DEC2 )/2. XMERI=(DZ1+DZ2)/2.D0 ULAT(N)=DECMIN +(XMIC+REF+FLEV+XMERI)/RHO ALAT(N)=ULAT(N) V(N)=STALAT-ULAT(N) NST1(N)=IST1(I) NST2(N)=IST2(I) CALL ARCS(DECMIN,IP,MP,SP) CALL ARCS(ULAT(N),IL,ML,SL) WRITE(6,2) N,NST1(N),NST2(N),IP,MP,SP,XMIC,FLEV,REF,IL,ML,SL 2 FORMAT(' ',I21,2I9,I4,I3,F7.3,F11.3,F12.3,F9.3,I6,I3,F7.3) C C REJECTING READING CONTAINING BLUNDERS C IF(DABS(V(N)).LT.REJ) GO TO 119 IREJ(N)=1 GO TO 100 119 IREJ(N)=2 DMN(N)=DM(I) 100 CONTINUE GO TO 93 101 SUM=0.0 K=0 DO 120 I=1,N IREQ=IREJ(I) GO TO (120,121),IREQ 121 SUM=SUM+ULAT(I) K=K+1 120 CONTINUE IF(K.GT.0) GO TO 122 XMEAN=STALAT STDS=0. STDM=0. GO TO 200 122 XMEAN1 =SUM/K C L=NUMBER OF REJECTED OBSERVATIONS. L=0 REJ =18.D0/RHO DO 123 I=1,N IREQ=IREJ(I) GO TO (123,124),IREQ 124 V(I)=ULAT(I)-XMEAN1 IF(DABS(V(I)).LT. REJ ) GO TO 123 IREJ(I)=1 L=L+1 123 CONTINUE IF(L.GT.0) GO TO 101 C K=NUMBER OF ACCEPTED OBSERVATIONS. K=0 REJ =3.D0/RHO DO 150 I=1,N IREQ=IREJ(I) GO TO (150,148),IREQ 148 V(I)=ULAT(I)-XMEAN1 IF(DABS(V(I)).LT. REJ ) GO TO 152 IREJ(I)=1 L=L+1 GO TO 150 152 K=K+1 DMM(K)=DMN(I) VM(K)=V(I) 150 CONTINUE IF(L.GT.0) GO TO 101 C C LEAST SQARES ADJUSTMENT TO OBTAIN THE MOST PROBABLE VALUE OF A TURN C OF THE MICROMETER AND THE MOST PROBABLE VALUE OF THE LATITUDE OF C THE STATION. C ATPA(1,1)=0. ATPA(1,2)=0. ATPA(2,1)=0. ATPA(2,2)=0. ATPW(1)=0. ATPW(2)=0. DO 114 J=1,K P(J)=1.D0 AMAT(J,1)= 1.0 AMAT(J,2)=-DMM(J) ATPA(1,1)=ATPA(1,1)+AMAT(J,1)*P(J)*AMAT(J,1) ATPA(1,2)=ATPA(1,2)+AMAT(J,1)*P(J)*AMAT(J,2) ATPA(2,1)=ATPA(1,2) ATPA(2,2)=ATPA(2,2)+AMAT(J,2)*P(J)*AMAT(J,2) ATPW(1)=ATPW(1)+VM(J)*P(J)*AMAT(J,1) ATPW(2)=ATPW(2)+VM(J)*P(J)*AMAT(J,2) 114 CONTINUE CALL JORDAN(ATPA,2,DET,2) CALL MATRIX(2,2,1,2,2,1,ATPA,ATPW,XMAT) CALL MATRIX(150,2,1,K,2,1,AMAT,XMAT,AX) C= XMAT(1) DR= XMAT(2) XMEAN=XMEAN1+C VV=0. M=0 DO 128 I=1,K VM(I)=AX(I)-VM(I) VV=VV+VM(I)*P(I)*VM(I)*RHO*RHO 128 CONTINUE NDF=K-2 EVF=VV/NDF ATPA(1,1)=ATPA(1,1)*EVF ATPA(2,1)=ATPA(2,1)*EVF ATPA(1,2)=ATPA(2,1) ATPA(2,2)=ATPA(2,2)*EVF STDM=DSQRT(ATPA(1,1)) DO 155 I=1,N IF(IREJ(I).EQ.1) M=M-1 IF(IREJ(I).EQ.1) GO TO 127 ALAT(I)=ULAT(I)+DMM(I+M)*DR GO TO 155 127 ALAT(I)=ULAT(I) 155 CONTINUE 200 WRITE(6,134) 134 FORMAT('1',36X,'VALUES OF ADJUSTED LATITUDES AND RESIDUALS',//,37X 1,'PAIR STAR1 STAR2 LATITUDE RES.',/,38X,'NO',30X,'(SECS)', 3/) L=N DO 138 I=1,N V(I)=ALAT(I)-XMEAN V(I)=V(I)*RHO 25 CALL ARCS(ALAT(I),ID(I),MI(I),SE(I)) WRITE(6,135) I,NST1(I),NST2(I),ID(I),MI(I),SE(I),V(I) 135 FORMAT(' ',36X,I3,2I7,I4,I3,F7.3,F8.3) IF(IREJ(I).NE.1) GO TO 138 L=L-1 WRITE(6,140) 140 FORMAT('+',80X,'REJECTED') 138 CONTINUE CALL ARCS(XMEAN,MDEG,MMIN,SECM) C=C*RHO DR=DR*RHO*2.D0 WRITE(6,27) MDEG,MMIN,SECM,L,STDM,C,DR,VV,NDF,EVF 27 FORMAT(' ',//,43X,'MEAN LATITUDE =',2I4,F8.3,5X,'FROM',I5,' PAIRS' 1,/,43X,'STANDARD DEV. OF MEAN = + ',F5.3,'SECS.',/,67X,'-',/,43X,' 2LATITUDE CORRECTION(AFTER ADJUSTMENT)=',F7.3,'SEC',/,43X,'CORRECTI 3ON TO MIC. VALUE=',F8.3,'SEC',/,43X,'VTPV=',F7.3,/,43X,'DEGREES OF 4 FREEDOM=',I3,/,43X,'ESTIMATED VAR. FACTOR=',F8.3,//) WRITE(6,156) 156 FORMAT(' ',////,30X,'VARIANCE-COVARIANCE MATRIX OF UNKNOWNS AFTER 1ADJUSTMENT',//) WRITE(6,154) ATPA 154 FORMAT(' ',47X,2F9.4) WRITE(6,92) 92 FORMAT('1') STOP END SUBROUTINE ARCS(R,ID,M,S) IMPLICIT REAL *8(A-H,O-Z) PI=3.14159265358979D0 DEGRAD=180./PI TWOPI=PI*2.D0 1 IF(R.LE.TWOPI) GO TO2 IR=R/TWOPI R=R-IR*TWOPI GO TO 1 2 DEG=R*DEGRAD ID=DEG FM=(DEG-ID)*60. M=FM S=(FM-M)*60. IF(DABS(S-60.).GT.1.D-12) GO TO 3 S=DABS(S-60.) M=M+1 IF(M.LT.60) GO TO 3 M=M-60 ID=ID+1 3 RETURN END SUBROUTINE MATRIX(N,K,M,NN,KK,MM,A,B,C) IMPLICIT REAL *8(A-H,O-Z) DIMENSION A(N,K),B(K,M),C(N,M) DO 1 I=1,NN DO 1 J=1,MM C(I,J)=0. DO 1 L=1,KK 1 C(I,J)=C(I,J)+A(I,L)*B(L,J) RETURN END C SUBROUTINE JORDAN(A,NORDER,DET,NDIM) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION A(NDIM,NDIM),Y(100),ICJ(100),IRR(100),JCC(100) IF(NDIM.LE.100) GOTO20 WRITE(6,10) 10 FORMAT(5X,'MATRIX LARGER THAN 100*100') RETURN 20 DET=1.D0 DO 205 K=1,NORDER PIVOT=0.D0 DO 100 I=1,NORDER DO 100 J=1,NORDER IF(K.EQ.1) GOTO50 KA=K-1 DO 30 IA=1,KA IF(I.EQ.IRR(IA)) GOTO100 IF(J.EQ.JCC(IA)) GOTO100 30 CONTINUE 50 IF(DABS(A(I,J)).LE.PIVOT) GOTO100 PIVOT=DABS(A(I,J)) IRR(K)=I JCC(K)=J 100 CONTINUE IF(PIVOT.GT.1D-51) GOTO120 WRITE(6,110) 110 FORMAT(5X,'PIVOT ELEMENT TOO SMALL') DET=0.D0 RETURN 120 IR=IRR(K) JC=JCC(K) PIVOT=A(IR,JC) DET=DET*PIVOT DO 130 J=1,NORDER 130 A(IR,J)=A(IR,J)/PIVOT A(IR,JC)=1.D0/PIVOT IF(NORDER.EQ.1) RETURN DO 205 I=1,NORDER AIJC=A(I,JC) IF(I.EQ.IR) GOTO205 A(I,JC)=-AIJC/PIVOT DO 200 J=1,NORDER IF(J.EQ.JC) GOTO200 A(I,J)=A(I,J)-AIJC*A(IR,J) 200 CONTINUE 205 CONTINUE DO 210 I=1,NORDER IR=IRR(I) JC=JCC(I) 210 ICJ(IR)=JC NA=NORDER-1 DO 220 I=1,NA IA=I+1 DO 220 J=IA,NORDER IF(ICJ(J).GE.ICJ(I)) GOTO220 TEMP=ICJ(J) ICJ(J)=ICJ(I) ICJ(I)=TEMP DET=-DET 220 CONTINUE DO 240 J=1,NORDER DO 230 I=1,NORDER IR=IRR(I) JC=JCC(I) 230 Y(JC)=A(IR,J) DO 240 I=1,NORDER 240 A(I,J)=Y(I) DO 260 I=1,NORDER DO 250 J=1,NORDER IR=IRR(J) JC=JCC(J) 250 Y(IR)=A(I,JC) DO 260 J=1,NORDER 260 A(I,J)=Y(J) RETURN END