C C ............................................................................ C C *** ASTRO *** C C VERSION JULY 1977 C C PROGRAMMED BY DEMITRIS DELIKARAOGLOU , UNB C FOR THE LRIS, SURVEY & MAPPING DIVISION, PEI C C A PROGRAM FOR THE DETERMINATION OF ASTRONOMIC AZIMUTHS ALONG WITH THEIR C ACCURACIES AS DERIVED FROM FIELD OBSERVATIONS ON POLARIS AT ANY HOUR ANGLE. C C THE PROGRAM UTILIZES THE UT OF OBSERVATION (TRANSFORMED FROM THE RECORDED C UTC TIME) TO COMPUTE THE HOUR ANGLE AND HENCE THE AZIMUTH OF POLARIS . C THE AZIMUTH OF THE LINE "STATION - RO" IS COMPUTED BY APPLYING THE CORRECTED C HORIZONTAL ANGLE BETWEEN THE STAR AND RO TO THE AZIMUTH OF POLARIS. C C IN ADDITION THE PROGRAM PERFORMS A CHI-SQUARE STATISTICAL TEST ON THE ESTIMA C TED VARIANCE , AND A NORMAL TEST OR ALTERNATIVELY A MAXIMUM-TAU TEST ON THE C COMPUTED RESIDUALS TO INDENTIFY POOR DATA SETS WHICH ARE TO BE FLAGGED FOR C REJECTION. C C .........................................................##................# C C INPUT DATA CARDS C C I - PROGRAM OPTION CARD CONTAINING OPTION CODES : C ALPHA = DESIRED SIGNIFICANCE LEVEL FOR THE STATISTICAL TESTS PERFORMED. C ILIST = LEVEL OF OUTPUT LISTING. C 1 - ONE LINE SUMMARY OF FINAL COMPUTED AZIMUTHS AND THEIR ACCURACIES. C 2 - FULL LISTING OF INPUT AND COMPUTED DATA. C 3 - SUMMARY OF COMPUTED INDIVIDUAL AZIMUTH VALUES AT EACH STATION C PLUS STATISTICAL TEST RESULTS. C 4 - COMPLETE INPUT AND OUTPUT LISTING. C ITEST = 1 FOR MAXIMUM-TAU TESTING OF THE RESIDUALS. C 2 FOR NORMAL TESTING ON THE RESIDUALS. C 3 FOR MAXIMUM-NORMAL TESTING ON THE RESIDUALS. C ACE = DESIGNED OR A PRIORI ASSUMED KNOWN VARIANCE FACTOR. C C OBSERVATION DATA CARD DECK AT A STATION : C C II - STATION HEADING CARD CONTAINING STATION PRINTOUT HEADING UP TO C 76 CHARACTERS ( FORMAT - 4X,76A1 ) C C III - OBSERVING STATION CARD CONTAINING : C IS1 , IS2 - OBSERVING STATION & RO NUMBERS ( FORMAT - I5 , 96 ) C ID,IM,IYR - DATE OF OBSERVATION (DAY,MONTH,YEAR) ( FORMAT - 293 , 95) C LATD, LATM,SLAT - LATITUDE OF OCCUPIED STATION (DEGREES,MINUTES, C SECONDS) POSITIVE NORTH ( FORMAT - I3,I2,F3.1,1X ) C LONGD,LONGM,SLONG - LONGITUDE OF OCCUPIED STATION (DEGREES,MINUTES, C SECONDS) POSITIVE EAST ( FORMAT - I3,I2,F3.1,1X) C FOR WEST LONGITUDES SET - LONGD ONLY C C IV - OBSERVATION DATA CARDS C C 1 - POLARIS CARD : C IRAD,IRAM,RAS - RIGHT ASCENSION (HOURS,MINUTES,SECONDS) C ( FORMAT - 2I2,F3.1,1X ) C IDECD,IDECM,DECS - DECLINATION (DEGREES,MINUTES,SECONDS) C ( FORMAT - 2I2,F3.1,1X ) C IUTD,IUTM - UT STARTING TIME (HOURS,MINUTES) C (E.G. COORDINATED UNIVERSAL TIME AS RECEIVED FROM RADIO C SIGNALS) C IR1D,IR1M,R1S,IR2D,IR2M,R2S - VALUES OF R (AS TABULATED IN THE STAR C ALMANAC) FOR THE 6 HOUR INTERVAL PRIOR AND AFTER C THE STARTING TIME UT ( FORMAT - 2I2 ) C TCR - TIME CORRECTION ( (CT-UTC) + DUT1 ) ( FORMAT - F5.1 ) C C 2 - OBSERVATION DATA : C RO - MEAN HORIZONTAL READINGS TO THE RO (DIRECT & REVERSE) C PLS - MEAN HORIZONTAL READINGS TO POLARIS (DIRECT & REVERSE) C TIME - MEAN TIME INTERVAL (DIRECT & REVERSE) FROM THE C STARTING TIME UT C C FOR ANOTHER SET OF OBSERVATIONS AT SAME STATION REPEAT III - IV OTHERWISE C USE A BLANK CARD TO TERMINATE OBSERVATION DATA DECK AT THIS STATION. C C REPEAT FROM II - IV FOR DATA DECK AT ANOTHER STATION. C C USE A CARD WITH A DELIMITER >0 IN COL 1-4 TO INDICATE THE END OF ALL DATA C CARDS FOR THIS RUN . C C PROGRAM MODE C PROCESS DATA STATION BY STATION AS FOLLOWS - C READ DATA SETS AT STATION C COMPUTE SINGLE AZIMUTH VALUES OF LINE "STATION - RO" FROM EACH SET C COMPUTE OVERALL MEAN AZIMUTH FROM ALL SETS C PERFORM STATISTICAL TESTING FOR THE ASSESSMENT OF ALL OBSERVATIONS. C C ............................................................................ C IMPLICIT REAL*8(A-H,O-Z) DIMENSION IHEAD(76),RO(6),PLS(6),TIME(4) DIMENSION AZI(100),V(100),IREJ(100) COMMON /ARRAYS/ IEDIT(100,2),EDIT(100,2) COMMON /PRMOPT/ ILIST,ITEST,ALPHA,ACE PI = DARCOS(-1.0D0) C C READ PROGRAM OPTION CARD. READ(5,1000) ALPHA,ILIST,ITEST,ACE C C SET DEFAULT VALUES IF( ALPHA.EQ.0.0 ) ALPHA = 0.01 IF( ILIST.EQ.0 ) ILIST = 3 IF( ITEST.EQ.0 ) ITEST = 3 IF( ACE.EQ.0.0 ) ACE = 33.3 C C START PROCESSING DATA AT STATION C M = 1 C C READ STATION HEADING CARD. C 22 READ( 5,1001 ) LN,(IHEAD(I),I = 1,76) C IF( LN.NE.0 ) GO TO 23 N = 1 C C READ OBSERVING STATION CARD C 11 READ(5,1002) IS1,IS2,ID,IM,IYR,LATD,LATM,SLAT,LONGD,LONGM,SLONG IF( IS1.EQ.0 ) GO TO 1 DLAT = LATD + ( LATM + SLAT/60. )/60. IF( LONGD.LT.0 ) DLONG = LONGD - ( LONGM + SLONG/60. )/60. C C READ ASTRO DATA CARDS C READ(5,1003) IRAD,IRAM,RAS,IDECD,IDECM,DECS,IUTD,IUTM, * IR1D,IR1M,R1S,IR2D,IR2M,R2S,TCR READ(5,1004) (RO(I),I=1,3),(PLS(I),I=1,3),(TIME(I),I=1,2), * (RO(I),I=4,6),(PLS(I),I=4,6),(TIME(I),I=3,4) C C C PLSD,PLSR - MEAN HORIZONTAL CIRCLE READINGS (DIRECT & REVERSE) TO POLARIS C PLSD = PLS(1)+( PLS(2) + PLS(3)/60. )/60. PLSR = PLS(4)+( PLS(5) + PLS(6)/60. )/60. C C C ROD , ROR - MEAN HORIZONTAL CIRCLE READINGS TO THE REFERENCE OBJECT C ROD = RO(1) + ( RO(2) + RO(3)/60. )/60. ROR = RO(4) + ( RO(5) + RO(6)/60. )/60. C C HCRD = PLSD - ROD HCRR = PLSR - ROR IF( HCRD.LT.0. ) HCRD = HCRD + 360. IF( HCRR.LT.0. ) HCRR = HCRR + 360. C C HCRM - (FINAL) REDUCED HORIZONTAL ANGLE BETWEEN POLARIS AND RO C HCRM = ( HCRD + HCRR )/2. C C C TIMED,TIMED - MEAN TIME INTERVALS (DIRECT & REVERSE) FROM STARTING TIME UT C TIMED = ( TIME(1) + TIME(2)/60. )/60. TIMER = ( TIME(3) + TIME(4)/60. )/60. C TIMEM = ( TIMED + TIMER )/2. C TCR = TCR/3600. C C C UT - UT OF OBSERVATION C UT = IUTD + IUTM/60. + TIMEM + TCR IF( UT.GT.24. ) UT = UT -24. C C C INTERPOLATE FOR R AT THE UT OBSERVATION TIME C R1 = IR1D + ( IR1M + R1S/60. )/60. R2= IR2D + ( IR2M + R2S/60. )/60. C DR21 = ( R2 - R1 )/6. DIF = DMOD( UT , 6.D0 ) C DR = DR21*DIF C AA=UT/6.D0 IAA = IDINT(AA) IAA = 6 * IAA C CALL DEG( HCRM,IH1,IH2,H3 ) CALL DEG( UT,JUTD,JUTM,UTS ) C C GST - GREENWICH SIDEREAL TIME C GST = UT + R1 + DR C C DR = DR*3600. DLONG = DLONG/15. C C DLST - LOCAL SIDEREAL TIME C DLST = GST + DLONG C CALL DEG( GST,IGD,IGM,GS ) CALL DEG( DLONG,LD,LM,SL ) CALL DEG( DLST,ILD,ILM,SLS ) IF( LD.GT.0 ) GO TO 10 LM = -LM SL = -SL C C RA, DECL - RIGHT ASCENSION AND DECLINATION OF POLARIS C 10 RA = IRAD + (IRAM + RAS/60.)/60. DECL = IDECD + (IDECM + DECS/60.)/60. DECL = DECL*PI/180. C C C H - HOUR ANGLE OF POLARIS C H = DLST - RA IF( H.LT.0.0 ) H = 24. + H CALL DEG( H,IHD,IHM,HS ) H = (H*15.)*PI/180. C DLAT = DLAT*PI/180. C C C AZPLS - AZIMUTH OF POLARIS C AZPLS = DATAN( DSIN(H)/(DSIN(DLAT)*DCOS(H) - * DCOS(DLAT)*DTAN(DECL) ) ) AZPLS = AZPLS*180./PI IF( AZPLS.LT.0. ) AZPLS = AZPLS + 360. CALL DEG( AZPLS,IAZD,IAZM,AZS ) C C C AZRO - AZIMUTH OF THE RO C AZRO = AZPLS - HCRM IF( AZRO.LT.0. ) AZRO = AZRO + 360. CALL DEG( AZRO,IAROD,IAROM,AROS ) C GO TO (80,85,80,85),ILIST 85 WRITE(6,2001) IHEAD WRITE ( 6,2002 ) LATD,LATM,SLAT,IS1,LONGD,LONGM,SLONG,IS2, * ID,IM,IYR,IRAD,IRAM,RAS,IDECD,IDECM,DECS,IUTD,IUTM WRITE(6,2003) N,(RO(I),I=1,6),(PLS(I),I=1,3),(TIME(I),I=1,2), * (PLS(I),I=4,6),(TIME(I),I=3,4) WRITE(6,2004) IH1,IH2,H3,JUTD,JUTM,UTS,IAA,IR1D,IR1M,R1S,DR WRITE(6,2005) IGD,IGM,GS,LD,LM,SL,ILD,ILM,SLS WRITE(6,2006) IRAD,IRAM,RAS,IHD,IHM,HS WRITE(6,2007) IAZD,IAZM,AZS,IH1,IH2,H3,IAROD,IAROM,AROS,IS1,IS2 C 80 AZI(N) = AZRO IEDIT(M,1) = IS1 IEDIT(M,2) = IS2 N = N + 1 C C TRANSFER TO STATEMENT 11 FOR NEXT DATA SET AT STATION. C GO TO 11 1 N = N - 1 IF( ILIST.EQ.2 ) GO TO 88 CALL RESIDU( M,N,AZI ) 88 M = M + 1 GO TO 22 23 M = M-1 GO TO (82,87,87,82),ILIST 82 CALL LIST(M) 87 STOP C C FORMAT STATEMENTS. C 1000 FORMAT(F4.2,2I4,F8.2) 1001 FORMAT( I4,76A1) 1002 FORMAT( I5,I6,2I3,I5,2(I3,I2,F3.1,1X) ) 1003 FORMAT( 2(2I2,F3.1,1X),2I2,2(1X,2I2,F3.1),F5.1 ) 1004 FORMAT ( 2( F3.0,F2.0,F3.1,F4.0,F2.0,F3.1,F3.0,F3.1,1X) ) 2001 FORMAT( '1',12X,80A1) 2002 FORMAT( '0',T40,'LATITUDE : ',2I3,F5.1,/,T10,'STATION : ',I5,/, * T40,'LONGITUDE : ',2I3,F5.1,//,T10,'R.O. : ',I5,//,T10, * 'DATE : ',2I3,I5,//,T40,'RIGHT ASCENSION :',2I3,F5.1,/,T10, * 'STAR : POLARIS',/,T40,'DECLINATION :',2I3,F5.1,//,T10, * 'START TIME :',2I3,' UT',/ ) 2003 FORMAT( '0',T40,'SET NO.',I3,//,T34,'MEAN READINGS',//,T12,'HCR', * T23,'DIRECT',T35,'TIME INT''VAL',5X,'REVERSE',5X,'TIME INT''VAL', * //,T12,'RO',6X,2F4.0,F5.1,16X,2F4.0,F5.1,//,T10,'POLARIS ', * 2F4.0,F5.1,3X,F4.0,F5.1,4X,2F4.0,F5.1,3X,F4.0,F5.1,/ ) 2004 FORMAT(//,T10,'MEAN REDUCED HCR : ',2I3,F5.1,///,T10, * 'UT OF OBSERVATION : ',2I3,F5.1,//,T10,'R AT ',I2,' HOURS UT', * T30,': ',2I3,F5.1,//,T10,'DELTA R',T30,': ',F5.1,///) 2005 FORMAT( T10,'GST',T30,': ',2I3,F5.1,//,T10,'LONGITUDE',T30,': ', * 2I3,F5.1,///,T10,'LST',T30,': ',2I3,F5.1,//) 2006 FORMAT(T10,'RA',T30,': ',2I3,F5.1,//,T10,'HOUR ANGLE OF STAR : ', * 2I3,F5.1,//) 2007 FORMAT(T10,'AZIMUTH OF STAR : ',2I3,F5.1,//,T10, * 'HORIZONTAL ANGLE : ',2I3,F5.1,///,T10,'*** TRUE AZIMUTH IS : *',2I3,F5.1,' FROM STN',I6,' TO STN',I6 ) END SUBROUTINE RESIDU( M,N,AZI ) C C PERFORMS A UNIVARIATE ASSESSMENT OF THE COMPUTED INDIVIDUAL AZIMUTH C VALUES AT STATION FOR QUALITY CONTROL C C INPUT PARAMETERS C C AZI - VECTOR OF INDIVIDUAL AZIMUTH VALUES COMPUTED FROM N DETERMINATIONS C N - LENGHT OF VECTOR AZI C M - ROW DIMENSION OF WORKING ARRAY FOR FINAL PRINTOUT C C SUBROUTINE MODE - C PUT TOGETHER INDIVIDUAL AZIMUTH VALUES TO FORM A "UNIVARIATE SET" C COMPUTE MEAN AND ESTIMATED VARIANCE C PERFORM A DISPERSION TEST ON THE VARIANCE C COMPUTE AND ANALYZE RESIDUALS C COME OUT WITH REPRESENTATIVE VALUE OF SET ( E.G. FINAL MEAN ) AND C ITS VARIANCE C IMPLICIT REAL*8( A-H,O-Z ) DIMENSION AZI(100),V(100),IREJ(100) COMMON/ARRAYS/IEDIT(100,2),EDIT(100,2) COMMON /PRMOPT/ ILIST,ITEST,ALPHA,ACE C C C COMPUTE MEAN AZIMUTH C AZM = 0.0 DO 20 I = 1,N 20 AZM = AZM + AZI(I) AZM = AZM / N C IS1 = IEDIT(M,1) IS2 = IEDIT(M,2) IF(ILIST.GE.3) WRITE (6,200) IS1,IS2,N C C COMPUTE AND ANALYSE RESIDUALS FOR QUALITY CONTROL C VV = 0.0 DO 21 I = 1,N V(I) = (AZM - AZI(I))*3600. V2 = V(I)**2 VV = VV + V2 IF( ILIST.LT.3 ) GO TO 21 CALL DEG(AZI(I),IA,MA,SA) WRITE(6,201) I,IA,MA,SA,V(I),V2 21 CONTINUE C SO = DSQRT(VV/DFLOAT(N-1)) STD = SO/DSQRT(DFLOAT(N)) SV = DSQRT( VV/DFLOAT(N) ) C IF( ILIST.LT.3 ) GO TO 29 CALL DEG(AZM,IAZ,MAZ,SAZ) WRITE(6,202) IAZ,MAZ,SAZ,SO,STD 29 NT = N NU = N-1 APVF = ACE 32 SIGMA = SO**2 C C PERFORM DISPERSION TEST ON THE VARIANCE C CALL ANALV( APVF,SIGMA,NU ) IF(ITEST.NE.1) GO TO 40 C C PERFORM TEST ON THE INDIVIDUAL RESIDUALS C CALL TAURE( NT,NU,ALPHA,TAU ) STDV = SV FCR = TAU GO TO 41 C 40 AA = 1. - ALPHA/2. IF(ITEST.EQ.3) AA = 1. - ( ALPHA/N )/2. SIG = DSQRT((DFLOAT(NU)/DFLOAT(N) )*APVF) C = XNORM( AA,0.0D0,1.0D0 ) BOJ$ STDV = SIG FCR = C C 41 GO TO (33,34,38), ITEST 33 WRITE(6,203) ALPHA GO TO 35 34 WRITE(6,213) ALPHA GO TO 35 38 WRITE(6,223) ALPHA 35 IRJ = 0 NN = N C C TEST ON THE INDIVUDUAL RESIDUALS C DO 22 I = 1,N IREJ(I) = 0 TESTV = DABS( V(I)/STDV ) - FCR IF( TESTV ) 22,5,5 5 IREJ(I) = I IRJ = IRJ + 1 AAZ = NN*AZM - AZI(I) NN = NN - 1 AZM = AAZ/NN 22 CONTINUE IF( IRJ.NE.0 ) GO TO 13 IF(ILIST.GE.3) WRITE(6,207) GO TO 12 C C IF REJECTIONS OCCURED RECOMPUTE MEAN AZIMUTH AND ACCURACY ESTIMATES C 13 IF(ILIST.GE.3) WRITE(6,208) VV = 0.0 DO 23 I = 1,N IF( I.NE.IREJ(I) ) GO TO 11 IF (ILIST.LT.3) GO TO 23 GO TO (36,37,39), ITEST 36 WRITE(6,204) I GO TO 23 37 WRITE(6,214) I GO TO 23 39 WRITE(6,224) I GO TO 23 11 VN = ( AZM - AZI(I) )*3600.0 V2 = VN*VN VV = VV + V2 IF(ILIST.LT.3) GO TO 23 CALL DEG(AZI(I),IA,MA,SA) WRITE(6,205) I,IA,MA,SA 23 CONTINUE SO = DSQRT(VV/DFLOAT(NN-1)) STD = SO/DSQRT(DFLOAT(NN)) SV = DSQRT( VV/DFLOAT(NN) ) IF(ILIST.LT.3) GO TO 12 CALL DEG(AZM,IAZ,MAZ,SAZ) WRITE(6,206) IAZ,MAZ,SAZ,SO,STD C C STORE FINAL MEAN AZIMUTH AND ITS STANDARD DEVIATION FOR SUMMARY PRINTOUT C 12 EDIT(M,1) = AZM EDIT(M,2) = STD C 200 FORMAT('1',25X,'SUMMARY',//,T6,' OBSERVED AZIMUTH OF LINE ', * I6,' - ',I6,/,T13,'FROM ',I2,' SETS OF DETERMINATIONS',//,5X, * 'SET NO.',6X,'OBS''D AZIMUTH',9X,'V',15X,'VV'/) 201 FORMAT(T8,I2,T19,I3,2X,I2,2X,F5.2,6X,F7.3, 9X,F7.3,/) 202 FORMAT(/,8X,'MEAN AZIMUTH',5X,'S.D. OF A SINGLE',6X,'S.D. OF THE $MEAN',/,9X,'OF THE LINE',4X,'AZIMUTH DETERMINATION',8X,'AZIMUTH', $ /,9X,2I3,F6.2,10X,F5.2,18X,F5.2,/) 203 FORMAT(////,35X,'MAXIMUM-TAU TESTING OF RESIDUALS',/, *35X,11('-'),1X,7('-'),1X,2('-'),1X,9('-'),/,20X,'REJECT SET IF STA *NTARDIZED RESIDUAL FAILS:',F5.2,' LEVEL OF TYPE I ERROR TEST',///) 213 FORMAT(////,38X,'NORMAL TESTING OF RESIDUALS',/, * 38X,6('-'),1X,7('-'),1X,2('-'),1X,9('-'),/,20X,'REJECT SET IF STA *NTARDIZED RESIDUAL FAILS:',F5.2,' LEVEL OF NORMAL TEST',///) 223 FORMAT(////,35X,'MAXIMUM-NORMAL TESTING OF RESIDUALS',/, *35X,14('-'),1X,7('-'),1X,2('-'),1X,9('-'),/,20X,'REJECT SET IF STA *NTARDIZED RESIDUAL FAILS:',F5.2,' LEVEL OF TYPE I ERROR TEST',///) 204 FORMAT(/,8X,I3,' TAU-TEST FAILED - SET REJECTED',/) 214 FORMAT(/,8X,I3,' NORMAL-TEST FAILED - SET REJECTED',/) 224 FORMAT(/,8X,I3,' MAXIMUM NORMAL-TEST FAILED - SET REJECTED',/) 205 FORMAT(8X,I3,7X,2I4,F7.2,/) 206 FORMAT(6X,'FINAL(MEAN) AZIMUTH',5X,'S.D. OF A SINGLE',7X,'S.D. OF *THE FINAL',/,9X,'OF THE LINE',8X,'AZIMUTH DETERMINATION',6X,'MEAN *AZIMUTH',//,9X,2I3,F6.2,10X,F5.2,18X,F5.2,/) 208 FORMAT(6X,'SETS KEPT OBS''D AZIMUTH',/) 207 FORMAT(37X,'NO REJECTIONS - ALL SETS KEPT',/) RETURN END SUBROUTINE ANALV( APVF,SIGMA,IDF ) C PERFORMS DISPERSION ( CHI-SQUARE ) TEST ON THE VARIANCE C INPUT PARAMETERS C APVF - A PRIORI VARIANCE FACTOR C SIGMA - ESTIMATED VARIANCE FACTOR C IDF - DEGREES OF FREEDOM IMPLICIT REAL *8(A-H,O-Z) COMMON /PRMOPT/ ILIST,ITEST,ALPHA IPR = ( (1.-ALPHA)*100. ) GO TO (15,15), ILIST WRITE(6,10) APVF,SIGMA,IDF 10 FORMAT( //,10X,'A PRIORI VARIANCE FACTOR 5',1X,'...',F6.2,//,10X,'ESTIMATED VARIANCE FACTOR',1X,'...',F8.4,/ 6/,10X,'DEGREES OF FREEDOM',1X,'...',I3,/////,33X,'CHI SQUARE TEST 7OF ESTIMATED VARIANCE FACTOR',/,33X,3('-'),1X,6('-'),1X,4('-'),1X, 82('-'),1X,9('-'),1X,8('-'),1X,6('-')) WRITE ( 6,12 ) ALPHA 12 FORMAT( /,34X,'GROUP OF OBSERVATIONS IS TO BE FLAGGED FOR',/,33X, * 'REJECTION IF FAILS :',F6.3,' LEVEL OF CHI SQUARE',/,37X,'TEST ON THE * ESTIMATED VARIANCE FACTOR') 15 AA = 1. - ALPHA/2. BB = ALPHA/2. AALPHA = CHISQ( AA,IDF ) BBETA = CHISQ( BB,IDF ) A=IDF*SIGMA/CHISQ( AA ,IDF) B=IDF*SIGMA/CHISQ( BB ,IDF) IF((IDF*SIGMA)/CHISQ( AA ,IDF) .LE. APVF) GO TO 4 5 CONTINUE IF(ILIST.GE.3) WRITE(6,70) IPR,APVF 70 FORMAT( ///,15X,'AT',I3,'% THE HYPOTHESIS THAT THE POPULATION VARI *ANCE EQUALS',F6.2,1X,'CAN BE REJECTED') RETURN 4 CONTINUE IF((IDF*SIGMA)/CHISQ( BB ,IDF) .LE. APVF) GO TO 5 IF(ILIST.GE.3) WRITE(6,71) IPR,APVF 71 FORMAT( ///,15X,'AT',I3,'% THE HYPOTHESIS THAT THE POPULATION VARI *ANCE EQUALS',F6.2,1X,'CAN NOT BE REJECTED') RETURN END SUBROUTINE TAURE( NT,NU,ALPH,CRTAU ) C COMPUTES THE REJECTION LEVEL FOR NORMALISED RESIDUALS FOR A GIVEN NUMBER OF C OBSERVATIONS , DEGREES OF FREEDOM AND DESIRED LEVEL OF TYPE I ERROR C PARAMETERS C NT - NUMBER OF OBSERVATIONS C NU - DEGREES OF FREEDOM C ALPH - DESIRED PROBABILITY OF TYPE I ERROR C CRTAU - CRITICAL VALUE ( MAX-TAU ) PRODUCED BY THE SUBROUTINE C REFERENCE : A.J. POPE (1976) - "THE STATISTICS OF RESIDUALS AND THE C DETECTION OF OUTLIERS" , U.S. DEPT. OF COMMERCE , NOAA TECHNICAL C REPORT NOS 65 NGS 1. IMPLICIT REAL*8(A-H,O-Z) DATA PI/ 3.1415926535898 / PD = 2. /PI S = 1. WNU = NU U = WNU -1. IF( U.EQ.0. ) GO TO 72 IF ( ALPH.EQ.0. ) GO TO 72 IF ( ALPH.LT.1. ) GO TO 10 CRTAU = 0. C RETURN C 10 Q = NT IF ( ALPH.GT.0.5 ) GO TO 19 AA = ALPH / Q DELT = AA DO 18 I = 1,100 XI = I DELT = DELT * ALPH * (( XI*Q - 1.)/(( XI+1.)*Q)) IF ( DELT.LT.1.D-14 ) GO TO 20 18 AA = AA + DELT 19 AA = 1. - (1.-ALPH)**(1./Q) 20 P = 1. - AA IF(U.EQ.1. ) GO TO 71 F = 1.3862943611199 - 2.*DLOG(AA) G = DSQRT(F) X = G - (2.515517 + 0.802853*G + 0.010328*F) $ / (1. + 1.432788*G + F*(0.189269 + 0.001308*G)) Y = X*X A = X*(1. + Y)/4. B = X*(3. + Y*(16. + 5.*Y))/96. C = X*(-15. + Y*(17. + Y*(19. + 3.*Y)))/384. E = X*(-945. + Y*(-1920. + Y*(1482. + Y*(776. + 79.*Y))))/92160. V = 1./U T = X + V*(A + V*(B + V*(C + E*V))) S = T/DSQRT(U + T*T) UM = U - 1. DELL = 1. DO 70 M = 1,50 H = 1. - S*S R = 0.0 IF ( DMOD(U,2.D0).EQ.0.0 ) GO TO 49 DD = DSQRT(H) N = 0.5*UM DO 45 I = 1,N Z = 2*I R = R + DD D = DD 45 DD = DD * H * (Z/(Z+1.)) R = PD*(R*S + DARSIN(S)) D = PD*D*UM GO TO 61 49 DD = 1. N = 0.5*U DO 55 I = 1,N Z = 2*I R = R + DD D = DD 55 DD = DD*H*((Z-1.)/Z) R = R*S D = D*UM 61 CONTINUE DEL = (P-R)/D IF( DABS( DEL/DELL ) .GT.0.5) GO TO 72 DELL = DEL S = S + DEL IF( DABS(DEL) .LT. 1.D-8 ) GO TO 72 70 CONTINUE GO TO 72 71 S =DSIN(P/PD) 72 CRTAU = S*DSQRT(WNU) RETURN END SUBROUTINE DEG(DECD,IDEG,MIN,SEC) C C CONVERTS DESIMAL DEGREES INTO DEGREES,MINUTES AND SECONDS. C IMPLICIT REAL *8(A-H,O-Z) C ID=IDINT(DECD) IDEG=ID DM=(DECD-ID)*60. IDM=IDINT(DM) MIN=IDM DS=(DM-IDM)*60. SEC=DS RETURN END SUBROUTINE LIST(M) C C PRINTS OUT ONE LINE SUMMARY LIST OF FINAL COMPUTED AZIMUTHS AND C THEIR ACCURACIES. C IMPLICIT REAL*8(A-H,O-Z) COMMON/ARRAYS/IEDIT(100,2),EDIT(100,2) C WRITE(6,101) 101 FORMAT('1',17X,'FROM',7X,'TO',7X,'OBS''D AZIMUTH',5X,'STANDARD',/, * 16X,'STATION STATION',23X,'DEVIATION') DO 10 I = 1,M CALL DEG ( EDIT(I,1),IAZ,MAZ,SAZ ) WRITE(6,100) IEDIT(I,1),IEDIT(I,2),IAZ,MAZ,SAZ,EDIT(I,2) 100 FORMAT('0',16X,I6,4X,I6,5X,I3,2X,I2,2X,F5.2,7X,F5.2) 10 CONTINUE RETURN END FUNCTION XNORM(ALF,XMEAN,SIG) C COMPUTES THE PERCENTILES OF THE NORMAL DISTRIBUTION N(XMEAM,SIG) C PARAMETERS C ALF - PROBABILITY INTEGRAL FROM NEGATIVE INFINITY TO XNORM C XMEAN - POPULATION MEAN C SIG - POPULATION STANDARD DEVIATION C XNORM - COPMUTED ABSCISSA VALUE OF N(XMEAN,SIG) CORRESPONDING TO C PROBABILITY ALF C ACCURACY BETTER THAN THAN 0.00045 0120 IMPLICIT REAL*8(A-H,O-Z) 0130 DIMENSION C(6) 0140 DATA C/2.515517D0,0.802853D0,0.010328D0,1.432788D0,0.189269D0, 0150 * 0.001308D0/ 0160 IF(ALF .GE. 1D0 .OR. ALF .LE. 0D0) GO TO 20 0170 SIGN = -1. 0180 P = ALF 0190 IF(ALF .LT. 0.5D0) GO TO 10 0200 SIGN = 1. 0210 P = 1. - ALF 0220 10 T = DSQRT(DLOG(1./P**2)) 0230 T2 = T * T 0240 T3 = T * T2 0250 XP = T - (C(1)+C(2)*T+C(3)*T2) / (1.+C(4)*T+C(5)*T2+C(6)*T3) 0260 XP = XP * SIGN 0270 XNORM = XMEAN + XP * SIG 0280 RETURN 0290 20 WRITE(6,1001) ALF 0300 1001 FORMAT(10X,'XNORM INPUT PROBABILITY=',E20.10) 0310 XNORM = 0. 0320 RETURN 0330 END 0340 FUNCTION CHISQ(A,N) C C COMPUTES THE PERCENTILES OF THE CHI-SQUARE DISTRIBUTION X(N) C C PARAMETERS C A - PROBABILITY INTEGRAL FROM ZERO TO CHISQ C N - NUMBER OF DEGREES OF FREEDOM C CHISQ - COMPUTED ABSCISSA VALUE OF X(N) CORRESPONDING TO PROBABILITY A C ACCURACY BETTER THAN 0.04 FOR N>1 IMPLICIT REAL*8(A-H,O-Z) CHISQ =N*(1-2./(9.*N)+XNORM(A,0D0,1D0)*DSQRT(2D0/(9.*N)))**3 RETURN END