FTN77,S $FILES 3,3 PROGRAM ANTPC(3,99) *********************************************************************** * * PROGRAM TO READ ANALYSE DATA COLLECTED BY TNFPC. * INFORMATION WRITTEN ON BY SUBROUTINE WINFO IN PROGRAM TNFPC * * This program writes data to final files * * THIS IS COPY OF NWNTS AS ON 1984/11/05 * ie can use other calibrators and get old nos from file * UPDATED <910118.0903> * * DOES NOT PRESERVE DVM READINGS. * *********************************************************************** * DOUBLE PRECISION DATE, STOBS * LOGICAL AOK, PROB, RPTFQ, HKEEP, ENDFIL, RPRINT, DUALBM, * FLXPNT, WVTB, CALSRS, WRTFNL, DONESX, DECDAT, DATF10, * FSTIM * INTEGER RESFIL(3), IDCB(144), RESLT(21), SORSNM(4), IPRAM(5), * ITIME(5), TMDCB(144), TMFILE(3), IBLNK(1), NAMFIN(3) * REAL RDFLX(11), ERFLX(11), PEAK(3), PEAKER(3), * DECOF(3), DECER(3) * COMMON / FLUXS / RDFLX, ERFLX, WAVE, SORSNM * DATA JREC / 0 /, JS / 0 /, LU / 1 /, * TMFILE / 2HTM, 2HFI, 2HLE / * * In/output unit is unit 1 CALL FSYSU (LU,LU) * CALL EXEC (11,ITIME,IYEAR) WRITE (LU,1010) ITIME(5),IYEAR 1010 FORMAT(//" PROGRAM ANTPC. Day ",I3,":",I5/ * 1X, 28"="/, * 3X,"Program to analyse pointing data."/) * WRITE (LU,1000) 1000 FORMAT(/" Program assumes that no more than 150 sets of "/ * " readings were taken."/ * " Check at end of TNFPC tracking program that"/ * " this number is not exceeded."/ * "If it is, increase all 200-dimensional arrays"/ * " in this program and recompile.") * 105 WRITE (LU,1100) 1100 FORMAT(/3X," Name of results file?_") READ (LU,1200) (RESFIL(L), L=1,3) 1200 FORMAT(3A2) * WRITE (LU,1300) (RESFIL(L), L=1,3) 1300 FORMAT(/3X," Reading results from disc file ",3A2, * /42X," OK (T)?_") * READ (LU,1500) AOK 1500 FORMAT(L6) IF ( .NOT. AOK ) GO TO 105 * C PRINT '(4X,"Date on file in decimal form (T or F) ? _")' C READ *, DECDAT DECDAT = .TRUE. C PRINT '(4X,"Data after D313/84 (T or F) ? _")' C READ *, DATF10 DATF10 = .TRUE. * C PRINT '(4X,"Data after D325/84 (T or F) ? _")' C READ *, FSTIM FSTIM = .TRUE. * C WRITE (LU,1350) C1350 FORMAT(/3X," Necessary to print out raw data ? Yes(T)_") C READ (LU,1500) RPRINT RPRINT = .FALSE. * C WRITE (LU,1360) C1360 FORMAT(/3X" Must NET FLUXES etc. be written out? Yes(T)_") C READ (LU,1500) FLXPNT FLXPNT = .FALSE. * WRITE (LU,1370) 1370 FORMAT(/3X," Writing results to final file ? Yes(T)_") READ (LU,1500) WRTFNL * * Open file, RESFIL. * 110 CALL OPEN (IDCB, IERR, RESFIL) IF ( IERR .GE. 0) GO TO 130 * * Error checking routine. * 120 CALL CLOSE (IDCB) WRITE (LU,1400) IERR 1400 FORMAT(" FMGR ERROR ",I4," in MAIN. CHECK STATUS. OK(T)? _") READ (LU,1500) AOK IF ( .NOT. AOK) GO TO 900 GO TO 110 * * Position pointer at begining of file. * 130 CALL POSNT (IDCB, IERR, 1, 1) IF ( IERR .LT. 0) GO TO 120 * * Read housekeeping, first resetting buffer. * 135 DO 133 IB = 1,21 RESLT(IB) = 2H 133 CONTINUE * CALL READF (IDCB, IERR, RESLT, 21, LEN) * JS = JS + 1 DONESX = .FALSE. * IF (( LEN .LE. 0 ) .OR. ( IERR .EQ. -012)) GO TO 150 IF (IERR .LT. 0) GO TO 120 GO TO 160 * 150 WRITE (LU,1600) 1600 FORMAT(/," END OF FILE ") GO TO 500 * * Assign to program parameters. * 160 CALL CODE READ (RESLT,1740) HKEEP 1740 FORMAT(L1) * * If line is line of housekeeping, continue. IF ( HKEEP ) GO TO 1750 * * Otherwise move back and write line of data. CALL POSNT (IDCB, IERR, -1) IF ( IERR .LT. 0 ) GO TO 120 GO TO 1850 * * If line is hkeeping, assign to program parameters * 1750 If (DECDAT) then If (DATF10) then If (FSTIM) then CALL CODE READ (RESLT,1706) HKEEP, NSORS, IBLNK, SORSNM, * DATE, KOUNT, CALSRS, SIDTIM else CALL CODE READ (RESLT,1705) HKEEP, NSORS, IBLNK, SORSNM, * DATE, KOUNT, CALSRS endif else CALL CODE READ (RESLT, 1700) HKEEP, NSORS,IBLNK,SORSNM, * DATE, KOUNT, CALSRS endif 1700 FORMAT(L1, I3, A2,4A2, F8.4, I2, L2) 1705 FORMAT(L1, I3, A2, 4A2, F10.6, I2, L2) 1706 FORMAT(L1, I3, A2, 4A2, F10.6, I2, L2, F12.6) * IYEAR = IDINT(DATE) IDAY = IDINT ((DATE - DBLE(IYEAR)) * 365.) IYEAR = IYEAR + 1900 * else CALL CODE READ (RESLT, 1710) HKEEP, NSORS, * IBLNK, SORSNM, IDAY, IYEAR, KOUNT, CALSRS 1710 FORMAT(L1, I3, A2, 4A2, I4, I4, I2, L2) DATE = FLOAT(IYEAR-1900) + FLOAT(IDAY)/365. * endif * C CALL NAMNO (NSORS, SORSNM, RAHR, *900,k) * * Find Hour Angle * DATSG = SNGL(DATE) STOBS = DBLE(SIDTIM) CALL HANGL (SORSNM, NSORS, DATE, HA, DATF10, FSTIM, * STOBS, *900,k) * * Want greater time resolution for 1144-379 (ie. >1 day) * If (.NOT. DECDAT) then If (NSORS .EQ. 123 .OR. NSORS .EQ. 916) then PRINT '("Source",I4,2X,4A2,". Enter UT of", * " observation in hh, mm : _")', * NSORS, SORSNM READ *, IHR, IMIN DTIM = (FLOAT(IHR) + FLOAT(IMIN)/60.)/(24.*365.) DATE = DATE + DTIM endif endif * * 1144-37 and 6 calibrators/controls common to both lists * IF (NSORS .EQ. 916) NSORS = 123 IF (NSORS .EQ. 944) NSORS = 144 IF (NSORS .EQ. 909) NSORS = 108 IF (NSORS .EQ. 915) NSORS = 119 IF (NSORS .EQ. 921) NSORS = 127 IF (NSORS .EQ. 929) NSORS = 132 IF (NSORS .EQ. 937) NSORS = 141 * * Write headers. IF ( .NOT. RPRINT ) GO TO 190 WRITE (LU,1900) 1900 FORMAT(//," No.",1X,"Source",3X,"Day:",6X,"No. freq.") * WRITE (LU,1800) NSORS, IBLNK, SORSNM, IDAY, IYEAR, KOUNT 1800 FORMAT(I3, A2, 4A2, I4, I5, I4) * * Write headers. 1850 IF ( .NOT. RPRINT ) GO TO 190 WRITE (LU,2100) 2100 FORMAT(/6X,"Freq",1X,"Place",7X,"Flux",19X,"DVM",/ * 6X," cm.",9X , 5X,"Jy ",12X,"Tube Source", * /6X,4"=",1X,5"=",1X,17"=",5X,"======= =======") * * Start of loop for each freq. * 190 DO 300 I = 1,KOUNT * * Start point for repeat of freq. * * READ DATA * 170 DO 310 J = 1,11 * DUALBM = ((WAVE .EQ. 3.4) .OR. (WAVE .EQ. 6.0)) * Only 5 readings taken at 13 and 18 single beam. IF (( J .LT. 6 ) .OR. DUALBM ) GO TO 320 * DO 325 JJ = 6,11 RDFLX(JJ) = 0.0 ERFLX(JJ) = 0.0 325 CONTINUE GO TO 330 * 320 CALL RESRD (WAVE, RDFLX(J),ERFLX(J),IDCB, RESLT, * LU, PROB, RPTFQ, HKEEP, ENDFIL, RESFIL, RPRINT, * WVTB, IFQ, JS, TUBE) * IF ( PROB ) GO TO 900 IF (ENDFIL) GO TO 500 * * If line was housekeeping, set read flux values to zero. IF ( HKEEP ) then RDFLX(J) = 0.0 ERFLX(J) = 0.0 WAVE = 0.0 * If reading not complete, ignore IF ((DUALBM .AND. (J .GT. 1 .AND. J .LT. 11)) .OR. * (.NOT. DUALBM .AND. (J .GT. 1 .AND. J .LT. 5))) * JS = JS - 1 GO TO 135 else endif * 310 CONTINUE * IF ( DONESX .AND. WAVE .EQ. 6.0) JS = JS + 1 * * Calculate net fluxes by subtracting off FN values. * 330 CALL NETTS (LU, FLXPNT, PEAK, PEAKER, AVPK, DECOF, JS, * RPTFQ, RAOFF, RAFER, PKMEAS, PKMSER, DECER, * DATSG, TUBE, IFQ) * If ( WAVE .EQ. 6.0 ) DONESX = .TRUE. * * Write info onto temporary file for later manipulation. * CALL TEMPF (LU, NSORS, SORSNM, WAVE, CALSRS, RPTFQ, * PEAK, PEAKER, AVPK, DECOF, DATSG, * PROB, TMDCB, TMFILE, JREC, JS, PKMEAS, * PKMSER, DECER, RAOFF, RAFER, HA, TUBE) IF (PROB) GO TO 900 * IF (RPTFQ) GO TO 170 * 300 CONTINUE * * Return for next source. GO TO 135 * 500 CALL CLOSE (IDCB) * * Correct for RA offset and calbrate. * CALL CALIB (LU, TMDCB, TMFILE, WRTFNL, RESFIL, * *900,k) * * Write data to mag tape * PRINT '(//4X,"Write data on disc file ",(3A2),"to ", * "mag tape ? Yes(T) _")', (RESFIL(L),L=1,3) READ *, AOK IF (AOK) CALL MAGTP (RESFIL, *900,k) * 900 CONTINUE * C WRITE (LU,9100) C9100 FORMAT(" Want to purge TEMPORARY FILE ? _") C READ (LU,1500) AOK C IF ( .NOT. AOK ) GO TO 910 * C CALL PURGE (TMDCB, IERR, TMFILE, 0, 22) C IF ( IERR .LT. 0 ) GO TO 120 * WRITE (LU,9999) 9999 FORMAT(///"*** END ANTPC ***"///) END * * SUBROUTINE TO READ RESULTS OFF DISC FILE RESFIL. * *********************************************************************** * SUBROUTINE RESRD (WAVE, SORSJY, ERROR, IDCB, RESLT, LU, * PROB, RPTFQ, HKEEP, ENDFIL, RESFIL, RPRINT, * WVTB, IFQ, JS, TUBE) * LOGICAL PROB, AOK, RPTFQ, HKEEP, ENDFIL, RPRINT, WVTB, CALCTB * INTEGER IDCB(1), RESLT(1), RESFIL(1), KPLACE(1) * PROB = .FALSE. ENDFIL = .FALSE. CALCTB = .FALSE. * * Read off file, first resetting buffer. * 100 DO 99 IB = 1,21 RESLT(IB) = 2H 99 CONTINUE * CALL READF (IDCB, IERR, RESLT, 21, LEN) IF (( LEN .LE. 0 ) .OR. (IERR .EQ. -012)) GO TO 110 IF (IERR .LT. 0) GO TO 120 GO TO 200 * * Error routine. 110 WRITE (LU,1100) 1100 FORMAT(" END OF FILE RESFIL ") ENDFIL = .TRUE. GO TO 900 * 120 CALL CLOSE (IDCB) WRITE (LU,1200) IERR 1200 FORMAT(" FMGR ERROR ",I4," in RESRD. CHECK STATUS. OK(T)? _") * READ (LU,1300) AOK 1300 FORMAT(L6) * IF ( AOK ) GO TO 180 * PROB = .TRUE. GO TO 900 * * File error corrected, so.... 180 CALL OPEN (IDCB, IERR, RESFIL) IF (IERR .LT. 0) GO TO 120 * * Assign read values to parameters. * 200 CALL CODE READ (RESLT,1430) HKEEP 1430 FORMAT(L1) * * If hkeeping, move back and leave subroutine * IF (HKEEP) GO TO 250 * CALL CODE READ (RESLT,1440) HKEEP, WVTB 1440 FORMAT(2L1) * IF ( .NOT. HKEEP ) GO TO 300 * 250 CALL POSNT (IDCB, IERR, -1) IF ( IERR .LT. 0 ) GO TO 120 GO TO 900 * * Read off wavelength and tube DVM values. * 300 IF ( .NOT. WVTB ) GO TO 310 C CALL POSNT (IDCB, IERR, -1) C IF (IERR .LT. 0) GO TO 120 C CALL READF (IDCB, IERR, RESLT, 21, LEN) CALL CODE READ (RESLT, 1450) HKEEP, WVTB, WAVE, AVTUB 1450 FORMAT(L1, L1, F6.1, G12.6) * CALCTB = .TRUE. * GO TO 100 * 310 CALL CODE READ (RESLT, 1400) HKEEP, WVTB, KPLACE, SORSJY, ERROR, * RPTFQ, AVDVM 1400 FORMAT(L1, L1, 1A2, G12.6, G12.6, L2, G12.6) * If (CALCTB) then If (WAVE .EQ. 3.4) IFQ = 1 If (WAVE .EQ. 6.0) IFQ = 2 If (WAVE .EQ. 13.) IFQ = 3 If (WAVE .EQ. 18.) IFQ = 4 TUBE = AVTUB - AVDVM endif * C GO TO 100 * IF ( .NOT. RPRINT ) GO TO 900 WRITE (LU,1500) WAVE, KPLACE, SORSJY, ERROR, AVTUB, AVDVM 1500 FORMAT(4X,F6.1,2X, 1A2,2X, G10.4,"+-", 3G10.4) * 900 RETURN END * * SUBROUTINE TO COMPUTE NET FLUXES AND OTHER PEARLS OF WISDOM. ********************************************************************** * SUBROUTINE NETTS (LU, FLXPNT, PEAK, PEAKER, AVPK, DECOF, JS, * RPTFQ, RAOFF, RAFER, PKAV, PKAVER, DECER, * DATE, TUBE, IFQ) * LOGICAL DUALBM, FLXPNT, RPTFQ * INTEGER TPLACE(7), SORSNM(4) * REAL AVFN(3), AVFNER(3), NETFLX(11), NETER(11), * RDFLX(11), ERFLX(11), PEAK(3), PEAKER(3), DECOF(3), * HPAV(2), HPAVER(2), HPRAT(2), * DECER(3), TUFLUX(4), ERFNDV(4), DE(4), JYDVM, * NETDVM, NTERDV * COMMON / FLUXS / RDFLX, ERFLX, WAVE, SORSNM * DATA TPLACE / 2HON, 2HHP, 2HHP, 2HON, 2HHP, 2HHP, 2HCE /, * TUFLUX / 175., 24., 35., 140.0 / * If (DATE .LT. 83.6466) then TUFLUX(1) = 167. TUFLUX(2) = 22.8 TUFLUX(3) = 32.8 endif * DUALBM = (( WAVE .EQ. 3.4) .OR. (WAVE .EQ. 6.0)) * * Average of first two FN's. AVFN(1) = (RDFLX(1) + RDFLX(5))/2. ERROR = 0.5 * SQRT( (ERFLX(1)*ERFLX(1)) * + (ERFLX(5)*ERFLX(5)) ) AVFNER(1) = ERROR * * Net flux of first peak. NETFLX(1) = RDFLX(3) - AVFN(1) ERROR = SQRT( (ERFLX(3)*ERFLX(3)) * + (AVFNER(1)*AVFNER(1)) ) NETER(1) = ERROR * * Net HP fluxes. NETFLX(2) = RDFLX(2) - AVFN(1) ERROR = SQRT( (ERFLX(2)*ERFLX(2)) * + (AVFNER(1)*AVFNER(1)) ) NETER(2) = ERROR * NETFLX(3) = RDFLX(4) - AVFN(1) ERROR = SQRT( (ERFLX(4)*ERFLX(4)) * + (AVFNER(1)*AVFNER(1)) ) NETER(3) = ERROR * * Ave. of centre FN's. AVFN(2) = (RDFLX(5) + RDFLX(7))/2. ERROR = 0.5 * SQRT( (ERFLX(5)*ERFLX(5)) * + (ERFLX(7)*ERFLX(7)) ) AVFNER(2) = ERROR * * Net centre flux NETFLX(7) = RDFLX(6) - AVFN(2) ERROR = SQRT( (ERFLX(6)*ERFLX(6)) * + (AVFNER(2)*AVFNER(2)) ) NETER(7) = ERROR * IF ( DUALBM ) GO TO 430 NETFLX(7) = 0.0 NETER(7) = 0.0 * * * Ave. of last FN's. 430 AVFN(3) = (RDFLX(7) + RDFLX(11))/2. ERROR = 0.5 * SQRT( (ERFLX(7)*ERFLX(7)) * + (ERFLX(11)*ERFLX(11)) ) AVFNER(3) = ERROR * * Net peak flux. NETFLX(4) = RDFLX(9) - AVFN(3) ERROR = SQRT( (ERFLX(9)*ERFLX(9)) * + (AVFNER(3)*AVFNER(3)) ) NETER(4) = ERROR * * Net HP's. NETFLX(5) = RDFLX(8) - AVFN(3) ERROR = SQRT( (ERFLX(8)*ERFLX(8)) * + (AVFNER(3)*AVFNER(3)) ) NETER(5) = ERROR * NETFLX(6) = RDFLX(10) - AVFN(3) ERROR = SQRT( (ERFLX(10)*ERFLX(10)) * + (AVFNER(3)*AVFNER(3)) ) NETER(6) = ERROR * IF ( .NOT. DUALBM) GO TO 400 * Find average of all FN's. TAVFN = (RDFLX(1) + RDFLX(5) + RDFLX(7) * + RDFLX(11))/4. SMFNSQ = RDFLX(1)*RDFLX(1) + RDFLX(5)* * RDFLX(5) + RDFLX(7)*RDFLX(7) * + RDFLX(11)*RDFLX(11) AVSQFN = SMFNSQ/4. * GO TO 420 * * Averages for single beam systems- 400 TAVFN = AVFN(1) SMFNSQ = RDFLX(1)*RDFLX(1) + RDFLX(5) * *RDFLX(5) AVSQFN = SMFNSQ/2. * 420 DEVFN = SQRT ( AVSQFN - TAVFN*TAVFN ) * * Calculate the correct errors from the old * JYDVM = TUFLUX(IFQ) / ABS(TUBE) * DO 425 IE = 1,4 If (IE .EQ. 1) IEE = 1 If (IE .EQ. 2) IEE = 5 If (IE .EQ. 3) IEE = 7 If (IE .EQ. 4) IEE = 11 ERFNDV(IE) = ERFLX(IEE) / (JYDVM * SQRT(3.)) CALL MEDIN (ERFNDV(IE), IE, EMED, DE) 425 CONTINUE * ERTUBE = (DE(2) + DE(3)) / 2. * DO 426 IN = 1,7 * If (.NOT. DUALBM .AND. IN .GE. 4) GO TO 426 NETDVM = NETFLX(IN) / JYDVM NTERDV = NETER(IN) / JYDVM * CDIFF = (NTERDV * NTERDV) + * (ERFNDV(1) * ERFNDV(1) * * (NETDVM / TUBE) * * (NETDVM / TUBE)) - * (3. * ERTUBE * ERTUBE / 2.) If (CDIFF .GT. 0.) then CNETER = SQRT(CDIFF) NETER(IN) = CNETER * JYDVM endif 426 CONTINUE * * Find average of peak fluxes. PKAV = ( NETFLX(1) - NETFLX(4))/2. PKAVER = 0.5 * SQRT( (NETER(1)*NETER(1)) * + (NETER(4)*NETER(4)) ) * IF ( .NOT. DUALBM) PKAV = PKAV * 2. * * Std. dev. of FN fluxes expressed as % of peak flux. PCNTDV = ( DEVFN/PKAV ) * 100.0 * * HP averages. DO 440 L1= 1,2 L2 = 3*L1 - 1 L3 = L2 + 1 L4 = L1 * L1 HPAV(L1) = ( NETFLX(L2) + NETFLX(L3))/2. ERROR = 0.5 * SQRT( (NETER(L2)*NETER(L2)) * + (NETER(L3)*NETER(L3)) ) HPAVER(L1) = ERROR * * Ratios. HPRAT(L1) = HPAV(L1) / NETFLX(L4) IF ( .NOT. DUALBM ) HPRAT(2) = 0.0 440 CONTINUE * CERAT = NETFLX(7) / PKAV * IF ( .NOT. FLXPNT ) GO TO 500 * * Write out calculated values. * WRITE (LU,2000) (SORSNM(L), L=1,4), WAVE 2000 FORMAT(//" Source ",4A2," at",F6.1, * "cm.",/ * 48"=") * WRITE (LU,2100) (TPLACE(L), L = 1,7) 2100 FORMAT(/," Place:",4X,A2, 6(8X,A2)) * WRITE (LU,2200) (NETFLX(L), L = 1,7) 2200 FORMAT("Net flux:", 7G10.4) * WRITE (LU,2300) (NETER(L), L = 1,7) 2300 FORMAT( 1X,"+-error",1X, 7G10.4) * WRITE (LU,2400) TAVFN, DEVFN, PCNTDV 2400 FORMAT(/" Average of FN values - std dev given" * " as % of PEAK FLUX ie. shows drift ",/ * 8X, G10.4,"+-",F10.4," ie. ", F10.4,"%") * WRITE (LU,2500) PKAV, PKAVER, HPAV(1), HPAVER(1), * HPRAT(1), HPAV(2), HPAVER(2), HPRAT(2), * CERAT 2500 FORMAT(/3X," Average peak flux (Jy): ",G10.4,"+-",G10.4/ * 6X," Half power averages: ", G10.4,"+-",G10.4, * " Ratio HP/Peak:", F10.4,/ * 28X, G10.4,"+-",G10.4,15X,F10.4,/ * 6X," Ratio CE/Peak : ",F10.4) * 500 CALL GSOFF (LU, NETFLX, WAVE, PEAK, PEAKER, DECOF, AVPK, JS, * RPTFQ, RAOFF, RAFER, NETER, FLXPNT, DECER, * AVPKER, DATE, IFQ) * * Set limit on PKAV. Assume values>100 are erroneous. * IF ( PKAV .LT. 0.0 ) PKAV = 0.0 IF ( ABS(PKAV) .GE. 100. ) PKAV = 99.99 IF ( ABS(PKAVER) .GE. 10.) PKAVER = 9.999 * RETURN END * ********************************************************************* * * SUBROUTINE TO CALCULATE OFFSET IN DEC USING 3 POINTS * ON ASSUMED GAUSSIAN. * ********************************************************************* * SUBROUTINE GSOFF (LU, NETFLX, WAVE, PEAK, PEAKER, DECOF, AVPK, * JS, RPTFQ, RAOFF, RAFER, NETER, FLXPNT, DECER, * AVPKER, DATE, IFQ) * LOGICAL DUALBM, RPTFQ, FLXPNT * REAL SIGDEC(4), PEAK(3), DECOF(3), NETFLX(1), NETER(1), * CNPK(3), BMSEP(4), PEAKER(3), SIGASS(4), DCFLIM(4), * NETFD, DECER(3) * DATA SIGDEC / .044, .082, .162, .247 /, * SIGASS / .048, .082, .162, .247 /, * DCFLIM / .055, .140, .185, .250 /, * SLPFAC / 10.9 /, FRCSFC / 0.2028 /, * BMSEP / .213, .242, 0.0, 0.0 /, * FAC1 / 0.693147181 / * If (DATE .GT. 85.595 .AND. IFQ .EQ. 3) then ! Feed off axis SIGDEC(3) = .158 DCFLIM(3) = .180 endif * DUALBM = (( IFQ .EQ. 1 ) .OR. ( IFQ .EQ. 2 )) * DO 110 LL = 1,2 L1 = LL*LL * IF ( (LL .EQ. 2) .AND. ( .NOT. DUALBM ) ) GO TO 110 * * Check for read fluxes being wrong sign due to low flux * or error in FN readings * RAT1 = NETFLX(L1)/NETFLX(L1+1) RAT2 = NETFLX(L1)/NETFLX(L1+2) * IF ((RAT1 .GE. 0.0) .AND. (RAT2 .GE. 0.0)) GO TO 120 * IF ( .NOT. FLXPNT ) GO TO 100 WRITE (LU,1000) 1000 FORMAT(" Incorrect sign. Check. PEAK set to 0.0 ") 100 PEAK(1) = 0.0 PEAK(2) = 0.0 PEAKER(1) = NETER(1) PEAKER(2) = NETER(4) DECOF(1) = 0.0 DECOF(2) = 0.0 DECER(1) = 0.0 DECER(2) = 0.0 C LL = LL + 1 GO TO 115 * 120 V = (NETFLX(L1) + NETFLX(L1+2)) / NETFLX(L1+1) * * Require SQRT(V+1) - 1 > 0 ie V>0 * IF ( V .LE. 0.0 ) GO TO 100 * X = SIGDEC(IFQ)*SIGDEC(IFQ)/(SIGASS(IFQ)*2.*FAC1) Y = (EXP(2.*FAC1*SIGASS(IFQ)*SIGASS(IFQ)/ * (SIGDEC(IFQ)*SIGDEC(IFQ))))/4. Z = (EXP(FAC1*SIGASS(IFQ)*SIGASS(IFQ)/ * (SIGDEC(IFQ)*SIGDEC(IFQ))))/2. SW = SQRT(Y+V) OFFST = X * ALOG(SW - Z) OFSER = 0.5 * (X/SW) / (SW - Z) * V1 = NETFLX(L1) + NETFLX(L1+2) V1ER = SQRT( (NETER(L1)*NETER(L1)) + (NETER(L1+2) * *NETER(L1+2)) ) VFRCER = ((V1ER*V1ER)/(V1*V1)) + ((NETER(L1+1)* * NETER(L1+1))/(NETFLX(L1+1)*NETFLX(L1+1))) VER = V * SQRT(VFRCER) C WRITE (LU,1010) V, VER C1010 FORMAT(" V ", F10.4,"+-", F10.6) * OFSER = OFSER * VER * DOS = OFFST / SIGDEC(IFQ) * * DOS limited by overflow error. * IF ( DOS .GT. 11. ) DOS = 11. * PEAK(LL) = NETFLX(L1) * EXP(FAC1 * DOS * DOS) * EFDD = EXP(FAC1*DOS*DOS) NETFD = NETFLX(L1) * EFDD * 2.* FAC1 * DOS / SIGDEC(IFQ) PERSQ = (EFDD*EFDD*NETER(L1)*NETER(L1)) + * (NETFD*NETFD*OFSER*OFSER) PEAKER(LL) = SQRT(PERSQ) * DECOF(LL) = OFFST DECER(LL) = OFSER * * Check that dec offset is less than limit for applicability * of gaussian fit in RA, else set flux to zero * If (ABS(DECOF(LL)) .GT. DCFLIM(IFQ)) GO TO 100 * GO TO (111,112) LL 111 IF ( PEAK(LL) .LT. 0.0 ) GO TO 100 GO TO 110 112 IF ( PEAK(LL) .GT. 0.0 ) GO TO 100 * 110 CONTINUE * * Single beam..... 115 IF ( DUALBM ) GO TO 130 AVPK = PEAK(1) AVPKER = PEAKER(1) DECOF(2) = 0.0 DECER(2) = 0.0 GO TO 140 * 130 AVPK = ( PEAK(1) - PEAK(2) )/2. AVPKER = 0.5 * SQRT( (PEAKER(1)*PEAKER(1)) + (PEAKER(2) * *PEAKER(2)) ) * 140 IF ( .NOT. FLXPNT ) GO TO 190 * WRITE (LU,1100) (DECOF(J), DECER(J), J=1,2), * ( PEAK(J),PEAKER(J), J=1,2), * AVPK, AVPKER 1100 FORMAT(" DEC Offsets : Beam A ", F10.6,"+-", F10.6,2X, * " Beam B ", F10.6,"+-",F10.6,/ * " Peak fluxes ", 7X, 2(3X, F10.4,"+-", F8.4)/ * " AVERAGE : ", F10.4,"+-", F10.6) * 190 IF ( IFQ .NE. 2 ) GO TO 900 * Calculate offset of CENTRE. * C IF ( .NOT. RPTFQ ) GO TO 200 * * If 2 readings at 6cm, choose which is to be used. * C CALL TWOSX (LU, RAOFF, NETFLX(7), AVPK) C GO TO 900 * 200 SLOPE = -1.0 * SLPFAC * AVPK * C ERSLOP = 10. * AVPKER RAOFF = NETFLX(7) / SLOPE * FRCNET = NETER(7) / NETFLX(7) FRCPK = AVPKER / AVPK SQOFER = FRCNET*FRCNET + FRCPK*FRCPK + FRCSFC*FRCSFC RAFER = ABS(RAOFF) * SQRT(SQOFER) * IF ( AVPK .NE. 0.0 ) GO TO 900 RAOFF = 0.0 RAFER = 0.0 * 900 RETURN END * ********************************************************************* * * SUBROUTINE TO WRITE INFO ONTO TEMPORARY FILE * ********************************************************************* * SUBROUTINE TEMPF (LU, NSORS, SORSNM, WAVE, CALSRS, RPTFQ, * PEAK, PEAKER, AVPK, DECOF, DATE, * PROB, TMDCB, TMFILE, JREC, JS, PKMEAS, * PKMSER, DECER, RAOFF, RAFER, HA, TUBE) * LOGICAL CALSRS, RPTFQ, AOK, PROB * INTEGER SORSNM(4), TMDCB(144), TMFILE(3), TMBUF(82), * JTMSIZ(4), LBUF(400) * REAL PEAK(3), PEAKER(3), DECOF(3), DECER(3) * DATA JTMSIZ / 0,24,0,0 /, JTMTYP / 3 /, JTMSEC / 0 /, * JTMCR / 52 /, JDCBS / 0 /, JNUM / 1 / * CALL LGBUF (LBUF,400) * * Set limit on DECOF for writing onto file. * DO 300 LL = 1,2 IF ( DECOF(LL) .GE. 10. ) DECOF(LL) = 9.99999 300 CONTINUE * * Open temp file. * 400 JREC = JREC + 1 405 CALL OPEN (TMDCB, IERR, TMFILE, 0, JTMSEC, JTMCR, JDCBS) IF ( IERR .EQ. -006 ) GO TO 410 IF ( IERR .GE. 0 ) GO TO 420 * * Error routine. * 415 CALL CLOSE (TMDCB) WRITE (LU, 4100) IERR 4100 FORMAT(/" FMGR ERROR ",I4," in TEMPF. Check. OK?_") READ (LU,4000) AOK 4000 FORMAT(L6) * IF ( AOK ) GO TO 405 PROB = .TRUE. GO TO 900 * * Create file if doesn't exist. 410 CALL CRETS (TMDCB, IERR, JNUM, TMFILE, JTMSIZ, JTMTYP, * JTMSEC, JTMCR, JDCBS) IF (IERR .LT. 0) GO TO 415 C GO TO 405 * * Position at start of file. 420 CALL POSNT (TMDCB, IERR, JREC, 1) IF (IERR .LT. 0) GO TO 415 * * Info to buffer. * CALL CODE WRITE (TMBUF, 4200) NSORS, (SORSNM(L), L=1,4), WAVE, * CALSRS, RPTFQ, (PEAK(L), L=1,2), * (PEAKER(L), L=1,2), AVPK, * (DECOF(L), L=1,2), (DECER(L), L=1,2), * DATE, JS, * PKMEAS, PKMSER, RAOFF, RAFER, HA, * TUBE 4200 FORMAT(I3, 4A2, F6.1, 2L1, 5G12.6, 4F7.5, F10.6,I3, F5.2,F5.3, * 2F7.5, F7.3, G12.6) * * Write onto file. CALL WRITF (TMDCB, IERR, TMBUF, 82) IF (IERR .LT. 0) GO TO 415 * * CALL CLOSE (TMDCB) * 900 RETURN END * ********************************************************************* * * SUBROUTINE TO CALCULATE SECOND OF TWO RA OFFSET VALUES * IF TWO MEASUREMENTS WERE MADE AT 6 cm. ON THE SAME SOURCE. * ********************************************************************* * SUBROUTINE TWOSX (LU, RAOFF, NETFLX, AVPK) * LOGICAL AOK * REAL NETFLX * SLOPE = -10.0 * AVPK OFST2 = NETFLX / SLOPE * 110 WRITE (LU,1100) RAOFF, OFST2 1100 FORMAT(" 6cm readings repeated."/ * " The 2 values of RA offset calculated are :", * 2(3X,F9.5),/" Which is wanted ? "/ * " Enter 1 if first value wanted, 2 for second"/ * " Enter 0 to get average _") * READ (LU,*) JCHC JCHC = JCHC + 1 * GO TO (111, 222, 333) JCHC * 111 PROVOF = ( RAOFF + OFST2 )/2. GO TO 400 222 PROVOF = RAOFF GO TO 400 333 PROVOF = OFST2 * 400 WRITE (LU,4100) PROVOF 4100 FORMAT(" Offset is : ", F9.5, " OK (T) ? _") * READ (LU,4200) AOK 4200 FORMAT(L6) * IF ( .NOT. AOK ) GO TO 110 * RAOFF = PROVOF * RETURN END * ********************************************************************* * * SUBROUTINE TO CORRECT FOR RA OFFSET AND CALIBRATE MEASURED FLUXES * Reads from TMFILE. * ********************************************************************* * SUBROUTINE CALIB (LU, TMDCB, TMFILE, WRTFNL, * RESFIL, k,*) * PARAMETER (NSCS = 400) * CHARACTER CALNAM(20)*8 * LOGICAL AOK, CALFND, CALSRS, RPTFQ, DUALBM, WRTFNL, NEWCAL, * SPCRAT * INTEGER TMDCB(144), TMFILE(3), TMBUF(82), SORSNM(4),NUMCL(20), * IBLNK(1), NAMFIN(3), NWCLNM(4), RESFIL(3), LBUF(400), * NDECF(4) * REAL PEAK(3), PEAKER(3), DECOF(3), SIGMA(4), * BMSEP(4), CNPK(2), CALFLX(4,20), CALER(4,20), FRAT(4), * FRATER(4), CNPKER(2), OFLIMT(4), * AVDECF(4,2), * DCDEV(4,2), WVLEN(4), * DECER(3), * RAFARY(NSCS), RARYR(NSCS), DCFARY(4,NSCS,2), * DCARYR(4,NSCS,2), DCF(NSCS), DCFR(NSCS) C EMA DCFARY,DCARYR * DATA WVLEN / 3.4, 6.0, 13.0, 18.0 /, * SIGMA / .048, .082, .162, .247 /, * BMSEP / .216, .254, 0.0, .250 /, * OFLIMT/ .072, .123, .123, .250 /, * FAC1 / 0.693147181 /, * CALFND / .TRUE. /, * NDUMY / 0 /, NCAL / 0 /, * IBLNK / 2H / * CALL LGBUF(LBUF,400) * NWCALN = 0 NRAF = 0 DO 501 NZ = 1,150 RAFARY(NZ) = 0. RARYR(NZ) = 0. DCF(NZ) = 0. DCFR(NZ) = 0. DO 500 NFQ = 1,4 NDECF(NFQ) = 0 DO 505 K = 1,2 DCFARY(NFQ,NZ,K) = 0. DCARYR(NFQ,NZ,K) = 0. 505 CONTINUE 500 CONTINUE 501 CONTINUE * * Open temp file. * 100 CALL OPEN (TMDCB, IERR, TMFILE) IF (IERR .GE. 0 ) GO TO 120 * * Error routine. * 110 CALL CLOSE (TMDCB) WRITE (LU,1100) IERR 1100 FORMAT(" FMGR ERROR ", I4," in CALIB. OK (T) ? _") READ (LU,1200) AOK 1200 FORMAT(L6) * IF ( .NOT. AOK ) GO TO 800 GO TO 100 * * Position at start of file. * 120 CALL POSNT (TMDCB, IERR, 1, 1) IF (IERR .LT. 0) GO TO 110 * * Write headers. IF ( .NOT. CALFND ) GO TO 125 * WRITE (LU,1300) 1300 FORMAT(//" CALIBRATORS AND CONTROLS ") WRITE (LU,1360) 1360 FORMAT(" No.",2X," Source ",2X," Wave ",3X," Flux (Jy)", * 6X," Dec offsets ",2X," RA off ",2X,"HA") GO TO 130 * 125 WRITE (LU,1330) 1330 FORMAT(//,3X,"FINAL VALUES",/ * 1X," No.",2X,"Source",4X,"Wave",4X,"Meas. Flux", * 4X,"Corr. Flux",3X,"Dec offsets",5X,"RA off", * 2X,"HA", * /1X,4"-",2X,6"-",4X,4"-",4X,10"-",4X,10"-",3X, * 11"-",5X, 6"-",2X,2"-") * * Read data. 130 CALL READF (TMDCB, IERR, TMBUF, 82, LEN) IF ((IERR .EQ. -012) .OR. (LEN .LE. 0)) GO TO 300 IF (IERR .LT. 0) GO TO 110 * * Assign values from buffer to program parameters. * CALL CODE READ (TMBUF, 1400) NSORS, SORSNM, WAVE, CALSRS, RPTFQ, * (PEAK(L), L=1,2), (PEAKER(L), L=1,2), * AVPK, (DECOF(L), L=1,2), * (DECER(L), L=1,2), DATE,JS, * PKMEAS, PKMSER, RAOFF, RAFER, HA, * TUBE 1400 FORMAT(I3, 4A2, F6.1, 2L1, 5G12.6, 4F7.5, F10.6,I3, F5.2,F5.3, * 2F7.5, F7.3, G12.6) * * Introduce IFQ for easier manipulation * IF ( WAVE .EQ. 3.4 ) IFQ = 1 IF ( WAVE .EQ. 6.0 ) IFQ = 2 IF ( WAVE .EQ. 13.0) IFQ = 3 IF ( WAVE .EQ. 18.0) IFQ = 4 DUALBM = ( (IFQ .EQ. 1) .OR. (IFQ .EQ. 2) ) * If ( CALFND ) then IF ( .NOT. CALSRS .AND. (NSORS .NE. NWCALN)) GO TO 130 * * Put calibrator values into array for later manipulation. If ( NSORS .NE. NDUMY .OR. (IFQ .EQ. KDNFQ)) then * NCAL = NCAL + 1 NDUMY = NSORS KDNFQ = IFQ endif endif * * Correct for RA offset using value obtained from 6 cm. * CALL RACOR (RAOFF, RAFER, PKMEAS, PKMSER, IFQ, DUALBM, * PEAK, PEAKER, FLXAV, FLXER, WVLEN, DATE) * * Multiply by calibration factor. IF ( CALFND ) GO TO 250 * FINFLX = FRAT(IFQ) * FLXAV FINER = FRAT(IFQ) * FLXER * SQFR1 = (FLXER*FLXER) / (FLXAV*FLXAV) SQFR2 = (FRATER(IFQ)*FRATER(IFQ)) / (FRAT(IFQ)*FRAT(IFQ)) * FINER = FINFLX * SQRT( SQFR1 + SQFR2 ) * * Set limit to FINER as written onto file in F5.3 IF (FINER .GT. 99.99) FINER = 99.99999 * * If measurement repeated choose which is wanted. 250 IF (RPTFQ) WRITE (LU,2222) 2222 FORMAT(" TEST. RPTFQ TRUE AFTER AVERAG CALC ") * * If searching for CALIBRATORS don't write onto file. * IF ( CALFND ) GO TO 290 * IF ( .NOT. WRTFNL ) GO TO 260 * Write results onto file * CALL WFINL (NSORS, SORSNM, WAVE, DATE, FINFLX, FINER, * DECOF, DECER, RAOFF, HA, RESFIL, *900,k) * C CALL WFIN2 (LU, NSORS, SORSNM, WAVE, DATE, C * FINFLX, FINER, PKMEAS, PKMSER, DECOF, RAOFF, C * NAMFIN) * 260 WRITE (LU,3000) NSORS, IBLNK, SORSNM, WAVE, PKMEAS, PKMSER, * FINFLX, FINER, (DECOF(L), L=1,2),RAOFF, * HA 3000 FORMAT(I4, A2, 4A2, F6.1, 2(2X,F5.2,"+-",F5.3),F7.3,2(2X,F7.3), * 2X,F4.1) * * Set limit at 2 Jy for use of values to calc ave offset * If (FINFLX .GE. 2.0) then * If (RAOFF .NE. 0.0 .AND. WAVE .EQ. 6.0) then NRAF = NRAF + 1 If (NRAF .GT. NSCS) then PRINT '(//"Too many sources. Change parameter NSCS ", * "in subroutine CALIB and recompile")' k = 1 GO TO 900 endif RAFARY(NRAF) = RAOFF RARYR(NRAF) = RAFER endif * If (DECOF(1) .NE. 0.0 .OR. DECOF(2) .NE. 0.0) then NDECF(IFQ) = NDECF(IFQ) + 1 If (NDECF(IFQ) .GT. NSCS) then PRINT '(//"Too many sources. Change parameter NSCS ", * "in subroutine CALIB and recompile")' k = 1 GO TO 900 endif DO 510 K = 1,2 DCFARY(IFQ,NDECF(IFQ),K) = DECOF(K) DCARYR(IFQ,NDECF(IFQ),K) = DECER(K) 510 CONTINUE endif endif * GO TO 130 * 290 CALFLX(IFQ,NCAL) = FLXAV CALER(IFQ, NCAL) = FLXER NUMCL(NCAL) = NSORS * Assign 295 to label NSTMT = 295 295 WRITE (CALNAM(NCAL),'(4A2)', iostat=ierr, err=820) * (SORSNM(LS), LS=1,4) WRITE (LU, 3100) NSORS, IBLNK, SORSNM, WAVE, FLXAV, FLXER, * (DECOF(L), L=1,2), RAOFF, HA 3100 FORMAT(I4, A2, 4A2, F6.1, 3X, F5.2,"+-", F6.3, 3(2X,F7.3), * 2X,F4.1) * * Return to read next line. GO TO 130 * 300 IF ( .NOT. CALFND ) GO TO 800 * * If there are no calibrators in list, choose another * If (NCAL .EQ. 0) then PRINT '(/"No calibrators in list."/"Want to specify", * " another source as a calibrator? Yes(T) _")' READ *, NEWCAL * If (NEWCAL) then PRINT '("Name of new calibrator ? _")' READ '(4A2)', (NWCLNM(LN), LN=1,4) CALL NAMNO (NWCALN, NWCLNM, RAHR,*900,k) NDUMY = 0 else PRINT '("Want to specify cal ratio? Yes (T)_")' READ *, SPCRAT If (SPCRAT) then 310 DO 315 LF = 1,4 PRINT '("Enter cal ratio and error", * " for",(F5.1)," cm : _")', * WVLEN(LF) READ *, FRAT(LF), FRATER(LF) 315 CONTINUE * DO 320 LF = 1,4 PRINT '(/"Cal ratio for",(F5.1), * " cm is,"(F7.4)," +-", * (F7.4))', * WVLEN(LF), FRAT(LF), * FRATER(LF) 320 CONTINUE * PRINT '(12X,"OK (T)? _")' READ *, AOK IF (.NOT. AOK ) GO TO 310 else DO 330 LN=1,4 FRAT(LN) = 1.0 FRATER(LN) = 0.0 330 CONTINUE endif CALFND = .FALSE. endif GO TO 120 else * * Choose CALIRATION RATIO. * CALL FCALB (LU, NCAL, CALFLX, CALER, FRAT, * FRATER, NUMCL, CALNAM, NEWCAL) * CALFND = .FALSE. * * Return to find final corrected and calibrated fluxes. GO TO 120 endif * 800 CONTINUE * * Calculate averages,std dev * CALL WMEAN (RAFARY, RARYR, NRAF, AVRAOF, RADEV, AV, ER, DEV) * DO 520 NFQ = 1,4 DO 525 K = 1,2 DO 526 NZ = 1,NDECF(NFQ) DCF(NZ) = DCFARY(NFQ,NZ,K) DCFR(NZ) = DCARYR(NFQ,NZ,K) 526 CONTINUE CALL WMEAN (DCF, DCFR, NDECF(NFQ), AVDECF(NFQ,K), * DCDEV(NFQ,K), AV, ER, DEV) 525 CONTINUE 520 CONTINUE * WRITE (LU,5100) AVRAOF, RADEV, NRAF, * (WVLEN(NFQ), (AVDECF(NFQ,K), * DCDEV(NFQ,K), K=1,2), NDECF(NFQ), NFQ = 1,4) 5100 FORMAT(//15X,"Ave RA offset :", 4X, F7.3,"+-",F8.3,3X,I4/ * 15X,"Ave DEC offsets at :"/ * 4(20X,F6.1,"cm : ",2(3X, F7.3,"+-", F8.3),3X,I4/)) * GO TO 900 * * Error routine 820 CALL ERMSG (ierr, NSTMT, *900,k) GO TO label * 900 RETURN k END * ********************************************************************* * * Subroutine to correct for RA offset * ********************************************************************* * SUBROUTINE RACOR (RAOFF, RAFER, PKMEAS, PKMSER, IFQ, DUALBM, * PEAK, PEAKER, FLXAV, FLXER, WVLEN, DATE) * LOGICAL DUALBM, OK6 * REAL SIGMA(4), BMSEP(4), CNPK(2), CNPKER(2), PEAK(3), * PEAKER(3), WVLEN(4), SEPER(4), RAFLIM(4) * * 3.4 6.0 13.0 18.0 DATA SIGMA / .046, .081, .162, .247 /, * BMSEP / .213, .242, 0.0, 0.0 /, * RAFLIM/ .054, .093, .175, .250 /, * SEPER / .002, .006, 0.0, 0.0 /, * FAC1 / 0.693147181 / * If (DATE .GT. 86.945) then ! New 3.6 cm feeds BMSEP(1) = .195 SEPER(1) = 0.0 endif * If (DATE .GT. 85.595 .AND. IFQ .EQ. 3) then !Feed off axis SIGMA(3) = .166 RAFLIM(3) = .170 endif * If (IFQ .EQ. 2) then If (PEAK(1) .EQ. 0.0 .OR. PEAK(2) .EQ. 0.0) then OK6 = .FALSE. else OK6 = .TRUE. endif endif * DUMMY = RAOFF / SIGMA(IFQ) * * Limit DUMMY for EXP IF ( ABS(DUMMY) .GT. 11. ) DUMMY = 11. A = EXP( -FAC1 * DUMMY * DUMMY) * DYER = RAFER / SIGMA(IFQ) EMDF2 = 2.* FAC1 * DUMMY * A AER = ABS (EMDF2 * DYER) C PRINT '(" A ", F10.4,"+-", F10.6)', A, AER * * * Single beam..... If ( .NOT. DUALBM ) then CNPK(1) = PEAK(1) / A CNPK(2) = 0.0 CNRSQ = (AER*AER)/(A*A) + (PEAKER(1)* * PEAKER(1))/(PEAK(1)*PEAK(1)) CNPKER(1) = CNPK(1) * SQRT(CNRSQ) CNPKER(2) = 0.0 FLXAV = CNPK(1) FLXER = CNPKER(1) * else DMINUS = (RAOFF - BMSEP(IFQ) - SEPER(IFQ)) / SIGMA(IFQ) * IF ( ABS(DMINUS) .GT. 11. ) DMINUS = 11. B = EXP( -FAC1 * DMINUS * DMINUS) * EMFD2 = 2.* FAC1 * DMINUS * B BER = ABS(EMFD2 * DYER) * EXPD = (RAOFF - SEPER(IFQ)) / SIGMA(IFQ) D = EXP( -FAC1 * EXPD * EXPD) DER = ABS(2. * FAC1 * EXPD * D * DYER) * C PRINT '(" B, D ",2(F10.4,"+-", F10.6))', B, BER, D, DER * DPLUS = ( RAOFF + BMSEP(IFQ)) / SIGMA(IFQ) * IF ( ABS(DPLUS) .GT. 11. ) DPLUS = 11. C = EXP( -FAC1 * DPLUS * DPLUS) EMFD2 = 2.* FAC1 * DPLUS * C CER = ABS(EMFD2 * DYER) * EXPE = (RAOFF + SEPER(IFQ)) / SIGMA(IFQ) E = EXP( -FAC1 * EXPE * EXPE) EER = ABS(2. * FAC1 * EXPE * E * DYER) C PRINT '(" C, E ",2(F10.4,"+-", F10.6))', C, CER, E, EER * DMB = D - B DMBER = SQRT( (DER*DER) + (BER*BER) ) CME = C - E CMEER = SQRT( (CER*CER) + (EER*EER) ) DBRTSQ = (DMBER*DMBER) / (DMB*DMB) CERTSQ = (CMEER*CMEER) / (CME*CME) * CNPK(1) = PEAK(1) / DMB PKRAT = PEAKER(1)/PEAK(1) CNPKER(1) = CNPK(1) * SQRT( DBRTSQ + (PKRAT*PKRAT)) * CNPK(2) = PEAK(2) / CME PKRAT = PEAKER(2)/PEAK(2) CNPKER(2) = CNPK(2) * SQRT( CERTSQ + (PKRAT*PKRAT)) * FLXAV = ( CNPK(1) + CNPK(2) ) / 2. FLXER = 0.5 * SQRT( (CNPKER(1)*CNPKER(1)) + (CNPKER(2)* * CNPKER(2)) ) * endif * If (.NOT. OK6) then DUM = RAFLIM(2) / (2. * SIGMA(IFQ)) A2 = EXP ( -FAC1 * DUM * DUM) FLXER = FLXER / A2 endif * If (ABS(RAOFF) .GT. RAFLIM(IFQ)) then FLXAV = 0.0 FLXER = 0.0 endif * IF ( FLXAV .LT. 0.0 ) FLXAV = 0.0 * RETURN END * ********************************************************************* * * SUBROUTINE TO WRITE FINAL DATA ONTO FILE * ********************************************************************* * SUBROUTINE WFINL (NSORS, SORSNM, WAVE, DATE, FINFLX, FINER, * DECOF, DECER, RAOFF, HA, RESFIL, k,*) * CHARACTER NAMFIN*14 * INTEGER SORSNM(4), RESFIL(3), FRSFIL(3) * REAL DECOF(2), FDECOF(2), DECER(2) * DATA NAMFIN(1:1) / 'C' /, NAMFIN(7:14) / '::27:4:4' / * NWAVE = IFIX(WAVE*10.) IF (WAVE .EQ. 13.) NWAVE = 13 IF (WAVE .EQ. 18.) NWAVE = 18 NSORS = NSORS + 100 WRITE (NAMFIN(2:3),'(I2)',iostat=ierr, err=820) NWAVE WRITE (NAMFIN(4:6),'(I3)',iostat=ierr,err=820) NSORS * C PRINT '("Name ",(A14))', NAMFIN * * File operations * * Open file,if it exists * Assign 110 to label NSTMT=110 110 OPEN (10, file=NAMFIN, iostat=ierr, err=810, status='OLD') * * Find EOF Assign 120 to label NSTMT=120 120 DO WHILE (ierr .ne. -1) READ (10, 1100, iostat=ierr, err=820, end=150) FDATE,FLUX, * ERROR, (FDECOF(L),L=1,2), FRAOFF, FHA, * (FRSFIL(L),L=1,3) end do Assign 150 to label NSTMT=150 150 BACKSPACE(10, iostat=ierr, err=820) * * Write info to file * Assign 160 to label NSTMT=160 160 WRITE (10, 1100, iostat=ierr, err=820) DATE, FINFLX, FINER, * (DECOF(L),L=1,2), RAOFF, HA, * (RESFIL(L),L=1,3) * * Write EOF Assign 170 to label NSTMT=170 170 ENDFILE (10, iostat=ierr, err=820) * Assign 180 to label NSTMT=180 180 CLOSE (10, iostat=ierr, err=820) * GO TO 900 * * Error handling * * Check if error due to file not existing. * 810 If (ierr .eq. 462) then PRINT '("Creating file ",(A14))', NAMFIN OPEN (10, file=NAMFIN, iostat=ierr, err=820, status='NEW') * Assign 815 to label NSTMT = 815 815 WRITE (10, '(4A2)',iostat=ierr, err=820)(SORSNM(L),L=1,4) GO TO 120 endif * * General error checking. * 820 Call ERMSG (ierr, NSTMT, *900, k) GO TO label * 1100 FORMAT (F10.6, F5.2, F5.3, 4F6.4, 3A2) * 900 RETURN k END * ********************************************************************* * * Subroutine to display errors * ********************************************************************* * SUBROUTINE ERMSG (ierr, NSTMT, k,*) * LOGICAL AOK * PRINT '("Error ",(I4)," at statement",(I4))', ierr, NSTMT PRINT '("OK(T) or end program(F) ? _")' READ '(L1)', aok If (.NOT. aok) k=1 RETURN k END * ********************************************************************* * * SUBROUTINE TO CALCULATE RATIO OF MEASURED FLUX TO KNOWN * CALIBRATOR OR CONTROL FLUX. * ********************************************************************* * SUBROUTINE FCALB (LU, NCAL, CALFLX, CALER, FRAT, FRATER, NUMCL, * CALNAM, NEWCAL) * CHARACTER*8 CALNAM(20) * LOGICAL AOK, NEWCAL * INTEGER NWANT(20), NUMCL(20) * REAL CALRAT(20), CRATER(20), CALFLX(4,20), CALER(4,20), * TRUFLX(4,6), FRAT(4), FRATER(4), VAL(20), ERR(20) * * WAVE 3.4 6.0 13.0 18.0 * 0915-11 DATA TRUFLX / 8.1, 13.5, 27.3, 0.0, * 1648+05 * 6.36, 11.33, 27.9, 0.0, * 2203-18 * 3.38, 4.02, 5.48, 0.0, * 0500+01 * 1.50, 1.95, 2.32, 0.0, * 1421-49 * 3.60, 5.54, 7.51, 0.0, * Specified during running * 0.0, 0.0, 0.0, 0.0 / * WRITE (LU,1100) NCAL, (NUMCL(L), L=1,NCAL) 1100 FORMAT(/I4," Calibrators used. Source nos :: ",20I4) * * Loop for each freq. * DO 120 IFQ = 1,4 * IF ( IFQ .EQ. 1 ) WAVE = 3.4 IF ( IFQ .EQ. 2 ) WAVE = 6.0 IF ( IFQ .EQ. 3 ) WAVE = 13.0 IF ( IFQ .EQ. 4 ) WAVE = 18.0 * WRITE (LU,1200) WAVE 1200 FORMAT(/" Calibrators at ", F5.1," cm."/ * 1X, 24"*"/ * /" No. Meas. flux True flux Ratio "/ * " ==== ========== ========= ===== ") * * Loop for each source. DO 130 MCAL = 1,NCAL * IF (CALNAM(MCAL) .EQ. '0915-11 ') KCAL = 1 IF (CALNAM(MCAL) .EQ. '1648+05 ') KCAL = 2 IF (CALNAM(MCAL) .EQ. '2203-18 ') KCAL = 3 IF (CALNAM(MCAL) .EQ. '0500+01 ') KCAL = 4 IF (CALNAM(MCAL) .EQ. '1421-49 ') KCAL = 5 IF (NEWCAL) KCAL = 6 * If ( CALFLX(IFQ,MCAL) .EQ. 0.0 ) then CALRAT(MCAL) = 0.0 CRATER(MCAL) = 0.0 GO TO 125 endif * If (NEWCAL .AND. MCAL .LT. 2) then PRINT '("Enter TRUE FLUX for ",(A8)," at ", * (F5.1)," cm. :: _")', CALNAM(MCAL), * WAVE READ *, TRUFLX(IFQ,KCAL) endif * * Calibration ratio for this source. * 115 CALRAT(MCAL) = TRUFLX(IFQ,KCAL) / CALFLX(IFQ,MCAL) CRATER(MCAL) = CALRAT(MCAL) * CALER(IFQ,MCAL) / CALFLX(IFQ,MCAL) * 125 WRITE (LU,1400) NUMCL(MCAL), CALFLX(IFQ,MCAL), * TRUFLX(IFQ,KCAL), CALRAT(MCAL), CRATER(MCAL) 1400 FORMAT(1X,I4, 3(3X,F10.4),"+-", F6.4) * 130 CONTINUE * * Set NWANT to zero GO TO 210 * 150 WRITE (LU,1500) 1500 FORMAT(/" Which are wanted ? "/ * 18X," Enter 11 if there are no readings at this freq.", * /8X, " Enter 1 for first, 2 for second, etc, ", / * 10X," More than one no. entered causes those to be" * " averaged : _") READ *, NWANT C READ '(20I2)' , NWANT IF ( NWANT(1) .NE. 11 ) GO TO 155 FRAT(IFQ) = 1.0 GO TO 120 * * Calc average 155 DO 160 L = 1,20 NUM = NWANT(L) IF ( NUM .EQ. 0 ) GO TO 170 VAL(L) = CALRAT(NUM) ERR(L) = CRATER(NUM) 160 CONTINUE * 170 CALL WMEAN (VAL, ERR, (L-1), FRAT(IFQ), FRATER(IFQ), AV,ER,DV) * WRITE (LU,1700) WAVE, FRAT(IFQ), FRATER(IFQ) 1700 FORMAT(" Calibrator ratio for ",F5.1," cm. measurements" * " is :", F10.4,"+-",F6.4/ * 8X, " OK(T) ? _") * READ (LU,2000) AOK 2000 FORMAT(L6) * IF ( AOK ) GO TO 120 210 DO 180 LL = 1,20 NWANT(LL) = 0 180 CONTINUE GO TO 150 * 120 CONTINUE * RETURN END * ********************************************************************* * * SUBROUTINE TO WRITE FINAL DATA ONTO DISC FILE. * ********************************************************************* * SUBROUTINE WFIN2 (LU, NSORS, SORSNM, WAVE, DATE, * FINFLX, FINER, PKMEAS, PKMSER, DECOF, * RAOFF, NAMFIN) * LOGICAL AOK * INTEGER SORSNM(4), NAMFIN(3), FNDCB(144), FNBUF(42), IBLNK(1) * REAL DECOF(3) * DATA IBLNK / 2H / * * Open file 100 CALL OPEN (FNDCB, IERR, NAMFIN) IF ( IERR .GE. 0 ) GO TO 120 IF ( IERR .EQ. -006) GO TO 777 * * Error routine. 110 CALL CLOSE (FNDCB) WRITE (LU,1100) IERR 1100 FORMAT(" FMGR ERROR ", I4," in WFINL. Check. OK(T)?_") READ (LU,1200) AOK 1200 FORMAT(L2) * IF ( .NOT. AOK ) GO TO 900 GO TO 100 * * CREATE FILE 777 CALL CREAT (FNDCB, IERR, NAMFIN, 128, 3, 0, 15) IF ( IERR .LT. 0 ) GO TO 110 * * Position at EOF. * 120 CALL POSNT (FNDCB, IERR, 32766) IF ( IERR .EQ. -012 ) GO TO 140 IF ( IERR .LT. 0 ) GO TO 110 * * Data to buffer. * 140 CALL CODE WRITE (FNBUF, 1400) NSORS, IBLNK,(SORSNM(L), L=1,4), WAVE, * DATE, FINFLX, FINER, PKMEAS, PKMSER, * (DECOF(J), J=1,2), RAOFF 1400 FORMAT(I4, A2, 4A2, F6.1, F12.6, 2G12.6, F5.2, F5.3, 3F6.4) * * Write onto file. * CALL WRITF (FNDCB, IERR, FNBUF, 42) IF ( IERR .LT. 0 ) GO TO 110 * * Write EOF. * CALL WRITF (FNDCB, IERR, FNBUF, -1) IF ( IERR .LT. 0 ) GO TO 110 * * Close file. * CALL CLOSE (FNDCB) * 900 RETURN END * * ********************************************************************* * * WRITE DISC DATA FILE NF... TO MAG TAPE * ********************************************************************* * SUBROUTINE MAGTP (RESFIL, k,*) * CHARACTER DATFIL*10, DBUFR*42, TPBUF*29 * LOGICAL AOK, EOLF * INTEGER ITIME(5), RESFIL(3) * DATA DATFIL (7:10) / '::20' /, NCR / 20 /, MTOUT / 8 / * WRITE (DATFIL(9:10),'(I2)', iostat=ierr, err=820) NCR WRITE (DATFIL(1:6),'(3A2)', iostat=ierr, err=820) * (RESFIL(L), L=1,3) * Get time * CALL EXEC (11, ITIME, IYEAR) * GO TO 101 * 100 PRINT '("Data on cartridge no :? _")' READ *, NCR WRITE (DATFIL(9:10),'(I2)', iostat=ierr, err=820) NCR * PRINT '("Magtape LU no : _")' READ *, MTOUT * 101 PRINT '("Writing data from ",(A10)," to magtape",(I3),/,4X, * " at",(2I5,4I3))',DATFIL,MTOUT,IYEAR,(ITIME(L),L=5,1,-1) * PRINT '(9X,"OK (T)? _")' READ *, AOK IF ( .NOT. AOK ) GO TO 100 * * * Start tape with EOF mark to be compatible with rest * PRINT '("Is this a new tape, Yes(T)? _")' READ *, AOK If ( AOK ) then ENDFILE (MTOUT, iostat=ierr, err=820) GO TO 121 else PRINT '("Tape at end of last file ? Yes(T) _")' READ *, EOLF endif * * If just written into previous file ie. at end of 2X EOF, write data * If (EOLF) GO TO 115 * * Data on tape terminated by double EOF mark. Find it * 111 Do while (ierr .ne. -1) Assign 110 to label nstmt=110 110 READ (MTOUT, iostat=ierr, err=820, end=112) TPBUF end do 112 Assign 113 to label nstmt=113 113 READ (MTOUT, iostat=ierr, err=820, end=115) TPBUF GO TO 111 * 115 Assign 116 to label NSTMT=116 116 BACKSPACE (MTOUT, iostat=ierr, err=820) * * Write name, date as 1st record * 121 Assign 120 to label NSTMT = 120 120 WRITE (MTOUT, 1100, iostat=ierr, err=820) DATFIL, IYEAR, * (ITIME(L), L=5,1,-1) * * Dump file to tape, line by line * Assign 130 to label NSTMT = 130 130 OPEN (20, file=DATFIL, iostat=ierr, err=820, status='OLD') * Do while (ierr .ne. -1) Assign 132 to label NSTMT=132 132 READ (20, 1300, iostat=ierr, err=820, end=140) DBUFR * Assign 135 to label NSTMT=135 135 WRITE (MTOUT, 1300, iostat=ierr, err=820) DBUFR end do * * Write double EOF * 140 Assign 145 to label NSTMT=145 145 ENDFILE (MTOUT, iostat=ierr, err=820) ENDFILE (MTOUT, iostat=ierr, err=820) * GO TO 850 * 820 CALL ERMSG (ierr, NSTMT, *900,k) GO TO label * 1100 FORMAT (A10, I4, 5I3) 1300 FORMAT (A42) * 850 CONTINUE PRINT '("Rewind tape ? Yes(T) _")' READ *, AOK If (AOK) then Assign 855 to label nstmt=855 855 REWIND (MTOUT, iostat=ierr, err=820) endif * 900 RETURN k END * *********************************************************************** * * SUBROUTINE TO GET CORRECT SOURCE NO. FOR FINAL FILE NAME * ie. return to old nos. * *********************************************************************** * SUBROUTINE NAMNO (NSORS, SORSNM, RAHR,k,*) * CHARACTER NAMSRS*8, SNAMFL*10, FILSNM*8, IBLNK*2 * LOGICAL TWOTRY * INTEGER SORSNM(4) * SNAMFL = 'PNTING::16' TWOTRY = .FALSE. * 100 NSRSOL = NSORS * Assign 110 to label NSTMT = 110 110 WRITE (NAMSRS, '(4A2)', iostat=ierr, err=820) * (SORSNM(LS), LS=1,4) * C PRINT '("Change no of source ",(A8)," obs no ",(I4))', C * NAMSRS, NSRSOL * Assign 120 to label NSTMT = 120 120 OPEN (30, file=SNAMFL, iostat=ierr, err=820, status='OLD') * Assign 130 to label NSTMT = 130 Do while (NAMSRS .NE. FILSNM) 130 If (SNAMFL .NE. 'PNTING::16') then READ (30, 1300, iostat=ierr, err=820, end=200) FILSNO, IBLNK, * FILSNM, RH, RM, RS * RAHR = RH + (RM + RS/60.)/60. ! RA in decimal hours * else READ (30,1400, iostat=ierr, err=820, end=200) FILSNO, IBLNK, * FILSNM, SRA, SDEC RAHR = SRA/15. endif * end do * C PRINT '("Source ",(A8)," Observing no and file no resp", C * (2I6))', FILSNM, NSRSOL, FILSNO * NSORS = FILSNO GO TO 900 * 200 If (.NOT. TWOTRY) then C PRINT '("No match for source ",(A8)," observing no. ",(I4), C * " in file ",(A10))', C * NAMSRS, NSRSOL, SNAMFL SNAMFL = 'POINTC::16' TWOTRY = .TRUE. GO TO 100 else PRINT '("No match for either file. Enter RA of source ", * (A8),"in hh, mm _")', NAMSRS READ *, IHS, IMS RAHR = FLOAT(IHS) + (FLOAT(IMS))/60. GO TO 900 endif * * Error routine 820 CALL ERMSG (ierr, NSTMT, *900,k) GO TO label * 900 CLOSE (30) * RETURN k * 1300 FORMAT (I3,A2,A8, F2.0, F3.0, F5.1) 1400 FORMAT (I3, A2, A8, 2F9.3) END * * ********************************************************************** * * SUBROUTINE TO CONVERT BETWEEN SAST, UT AND ST * ********************************************************************** * SUBROUTINE TMCON(LU,ITYPE,TIMIN,HMIN,ID,IYR,KTPOUT,TIMOUT,HMOUT) * DOUBLE PRECISION TIMIN, TIMOUT, BASEOF, RATE, OFFJD, SLONG, * FAC1, FAC2, FAC3, FAC4, RJDNT0, SASTIM, UTIME, * STIME, GAST, DUMY, DECML, CHHMM * LOGICAL AOK, LEAPYR * INTEGER HMIN, HMOUT * C REAL JDNT0 * DATA BASEOF / 84 806.882 /, LSTDAY / 0 /, * RATE / 0.002 737 909 3 /, * OFFJD / 7051.0 /, SLONG / -1.8456667 /, * FAC1 / 6.67170278/, FAC2 / 0.0657098232 /, * FAC3 / 2433282.5 /, FAC4 / 2445334.5 / * * ITYPE is form of time coming IN. SAST ..... 1 * UT ..... 2 * ST ..... 3 * KTPOUT is reqd form of OUTPUT time. * Values as for ITYPE * * HMIN,HMOUT are set to 1 for input of form hh.mm, * 0 in decimal hrs. * * Calculations done in decimal form * IF ( HMIN .EQ. 1 ) TIMIN = DECML(TIMIN) * * Find Julian day no.; modified * LEAPYR = ( MOD(FLOAT(IYR)/4.,1.) .EQ. 0.0 ) If (LEAPYR) then LPY = 1 else LPY = 2 endif * If (ID .LE. (61-LPY)) then MNB = 0 else MNB = 1 endif * NTYR = IYR - 1960 * RJDNT0 = FLOAT( 367*NTYR - IFIX(7.* FLOAT(IYR+MNB)/4.) * + ID + LPY*MNB + 30) + OFFJD * * Do calcs depending on input time type. * GO TO (100,200,300) ITYPE * 100 SASTIM = TIMIN CALL LIM24(SASTIM, TMULT) * UTIME = SASTIM - 2. CALL LIM24(UTIME, TMULT) * GO TO 210 * * UT input * 200 UTIME = TIMIN CALL LIM24(UTIME, TMULT) * SASTIM = UTIME + 2. CALL LIM24(SASTIM, TMULT) * 210 GAST = FAC1 + FAC2*RJDNT0 + (1.+RATE)*UTIME STIME = GAST - SLONG CALL LIM24(STIME, TMULT) * GO TO 400 * * If ST entered * 300 STIME = TIMIN CALL LIM24(STIME, TMULT) * * Green app. ST * 310 GAST = STIME + SLONG * C WRITE (LU,3200) SLONG, GAST C3200 FORMAT(" Gren app ST from longitude ",F10.7,"is :",F10.4) * DUMY = FAC2*RJDNT0 UTIME = (GAST - FAC1 - DUMY) / (1. + RATE) * If (ABS(UTIME) .GT. 24.) then CALL LIM24 (UTIME, TMULT) STIME = STIME + 24.*TMULT GO TO 310 endif CALL LIM24(STIME, TMULT) * 330 SASTIM = UTIME + 2. CALL LIM24(SASTIM, TMULT) * 400 CONTINUE * * Write * C WRITE (LU,4100) SASTIM, UTIME, STIME C4100 FORMAT(/5X,"SAST", 6X,"UT", 6X,"ST",/ C * 3(3X,F6.2)) * * * Output time form as reqd. Set by value of KTPOUT * GO TO (511,522,533) KTPOUT * 511 TIMOUT = SASTIM GO TO 900 * 522 TIMOUT = UTIME GO TO 900 * 533 TIMOUT = STIME * 900 IF ( HMOUT .EQ. 0 ) GO TO 910 * * Else convert to hh.mm * ATIM = TIMOUT TIMOUT = CHHMM(ATIM) * 910 RETURN END * ********************************************************************** * * SUBROUTINE TO GET 0 < TIME < 24 * ********************************************************************** * SUBROUTINE LIM24 (XTIME, TMULT) * DOUBLE PRECISION XTIME * TMULT = 0.0 * 110 IF ( XTIME .LT. 24. ) GO TO 120 XTIME = XTIME - 24. TMULT = TMULT + 1. GO TO 110 * 120 IF ( XTIME .GT. 0.0 ) GO TO 130 XTIME = XTIME + 24. TMULT = TMULT + 1. GO TO 120 130 RETURN END * * *********************************************************************** * * Subroutine to calculate HOUR ANGLE from decimal date * *********************************************************************** * SUBROUTINE HANGL (SORSNM, NSORS, DATE, HA, DATF10, FSTIM, * STOBS, k,*) * DOUBLE PRECISION UTOBS, STOBS, UTCALC, DATE, DATNW * LOGICAL DATF10, FSTIM * INTEGER SORSNM(4) * * Date into year,day,UT * IYEAR = IDINT(DATE) DAY = 365. * (DATE - FLOAT(IYEAR)) IDAY = IFIX(DAY) IYEAR = IYEAR + 1900 UTOBS = 24. * (DAY - FLOAT(IDAY)) * C If (DATF10) then * * Convert UT to ST * C CALL TMCON (LU, 2, UTOBS, 0, IDAY, IYEAR, C * 3, STOBS, 0) C CALL LIM24 (STOBS,TMULT) * C else * If (.NOT. FSTIM) then PRINT '("Enter ST of obs of ",I4,2X,4A2," _")', NSORS, * (SORSNM(L),L=1,4) READ *, SHR, SMIN * STOBS = SHR + SMIN/60. endif * * Calculate UT from ST to get better precision for UT * CALL TMCON (LU, 3, STOBS, 0, IDAY, IYEAR, * 2, UTCALC, 0) CALL LIM24 (UTCALC,TMULT) * If (UTOBS .LT. 2. .AND. UTCALC .GT. 22.) IDAY = IDAY - 1 * * New decimal year from this UT DATNW = (UTCALC/24. + DBLE(IDAY))/365. + DBLE(IYEAR-1900) DATE = DATNW * C endif * * Source coords * CALL NAMNO (NSORS, SORSNM, RAHR, *900,k) * * HA = STOBS - RAHR IF ( HA .GT. 12. ) HA = HA - 24. IF ( HA .LT. -12.) HA = HA + 24. * HA = HA * 15. ! HA in degrees * C PRINT '("UTOBS, UTCALC, ST, RA, HA",5(2X,F8.6),/, C * "Filed date , New date",2(2X,F12.6))', C * UTOBS, UTCALC, STOBS, RAHR, HA, DATE, DATNW 900 RETURN k END ********************************************************************* * * FUNCTION TO GIVE DECIMAL TIME * ********************************************************************* * FUNCTION DECML (ATIM) * DOUBLE PRECISION ATIM, DECML * NTIM = IDINT(ATIM) ATIM = NTIM + ((ATIM - NTIM) * 100./60.) DECML = ATIM * RETURN END * ********************************************************************* * * FUNCTION TO CONVERT FROM DECIMAL TO HH.MM * ********************************************************************* * FUNCTION CHHMM (ATIM) * DOUBLE PRECISION ATIM, CHHMM * NTIM = IDINT(ATIM) ATIM = NTIM + ((ATIM - NTIM) * 60./100.) CHHMM = ATIM * RETURN END * *********************************************************************** * * Subroutine to calculate the weighted mean from arrays of * value and error. Also gives arithmetic mean. * *********************************************************************** * SUBROUTINE WMEAN (VALUE, ERROR, NDATA, WTMEAN, WTERR, AVERAG, * ARERR, DEVN) * PARAMETER (MXRD=150) * REAL VALUE(MXRD), ERROR(MXRD), WEIGHT(MXRD), WTVAL(MXRD) * IF (NDATA .EQ. 0) GO TO 900 * * Set sum values to zero * SUMWT = 0. SUMWV = 0. * SUMVAL = 0. SUMSQS = 0. SMSQER = 0. * DO 200 J = 1,NDATA * WEIGHT(J) = 1./(ERROR(J)*ERROR(J)) WTVAL(J) = WEIGHT(J) * VALUE(J) * SUMWV = SUMWV + WTVAL(J) SUMWT = SUMWT + WEIGHT(J) * SUMVAL = SUMVAL + VALUE(J) SUMSQS = SUMSQS + (VALUE(J) * VALUE(J)) SMSQER = SMSQER + (ERROR(J) * ERROR(J)) * 200 CONTINUE * WTMEAN = SUMWV/SUMWT WTERR = SQRT(1./SUMWT) * RNDATA = NDATA AVERAG = SUMVAL / RNDATA AVGSQD = AVERAG * AVERAG AVSMSQ = SUMSQS / RNDATA DIFF = AVSMSQ - AVGSQD If (DIFF .GE. 0.) then DEVN = SQRT (DIFF) else DEVN = 0.0 endif DEVN = DEVN / SQRT(RNDATA) ARERR = (SQRT(SMSQER)) / RNDATA * 900 RETURN END * ********************************************************************* * BLOCK DATA * INTEGER SORSNM * COMMON / FLUXS / RDFLX(11), ERFLX(11), WAVE, * SORSNM(4) * END * ********************************************************************* * * SUBROUTINE TO ARRANGE VALUES IN ASCENDING ORDER AND * CALCULATE THE MEDIAN. * Uses REAL variables. * ********************************************************************* * SUBROUTINE MEDIN (RDVM, ICOUNT, MEDVM, DIGVM) * REAL DIGVM(11), MEDVM * DO 110 J = 1,ICOUNT * * Rise thru' array until position of value is found. IF ( RDVM .LE. DIGVM(J) ) GO TO 130 * 110 CONTINUE * DIGVM(ICOUNT) = RDVM GO TO 900 * * Shift all subsequent values up the array. 130 DO 140 K = ICOUNT,J,-1 DIGVM(K+1) = DIGVM(K) 140 CONTINUE * DIGVM(J) = RDVM * 900 CONTINUE * MIDL = ICOUNT/2 MEDVM = DIGVM(MIDL) * RETURN END END$