Program FRITZ C Version 8.02 C C by David Mozurkewich 12 Sept 1991 C C This program reads the laser positions from the Mt Wilson C Interferometer data tapes with the snazzy new format. C It combines the old programs STAGE, AVER, COMPACT, VIS C Allows multiple baselines and any catalog you want. C C It outputs one color astrometry and visibility data. C seeing estimate is from laser data only. C C MODIFICATIONS: C laser rms is output instead of phi rms for seeing C calibration. C format version number is in vis file header. C c AQ and CAH: array parameter nlimit3, nlimit4 c limits (contains star numbers) c if limits(i+1)=9999, ask for new nlimit3 and -4 for each scan c for star limits(i). c IMPLICIT NONE REAL *8 CONVL, PI, RAD, LON, LAT C CONVL converts the LASER measurement, which is recorded C in double fringe units PARAMETER ( CONVL = 0.6329914D0 / 64.0D0 ) PARAMETER ( PI = 3.14159265358979D0 ) PARAMETER ( RAD = 180.D0 / PI ) PARAMETER ( LON = 118.059167D0 ) PARAMETER ( LAT = 34.2166944D0) INTEGER *4 MAXOBS, MAXCNT, MAXLOCK PARAMETER ( MAXOBS = 5000) PARAMETER ( MAXCNT = 600) PARAMETER ( MAXLOCK= 4096) CHARACTER *12 VERSN, VERSD PARAMETER( VERSN = 'VERSION 8.02' ) PARAMETER( VERSD = '1991 Sept 12' ) LOGICAL MOREDATA, NEWSTAR, LOSTLOCK, FAST, ALL INTEGER *4 IOBS, J, NBODY, TDTOFF INTEGER *4 K, ICNT, ICNT8, MINAVG INTEGER *4 STRTYR, STRTMON, STRTDAY, SPECT INTEGER *4 STRTSTAR, STRTSEC, STOPSEC, I, TEMP REAL *4 PMRA4, PMDEC4, PAR4, RV4, MAG4, BV, FDARK(4) REAL *4 TCOH, ANGDIA, TFIRST, TLAST REAL *8 LASER0 REAL *8 T0, TSUM, LSUM, LAVG, DJUL, UT1, GAST, TJD REAL *8 TAVG, RAA, DECA, DELAYC REAL *8 XSUM, YSUM, YYSUM, YAVG, RAM, DECM, MAG REAL *8 PMRA, PMDEC, PAR, RV, NPTS REAL *8 INTTIM, ETIME, TIME0 REAL *8 DELT, DLAS, LRMS, SLOPE, CONST C C Visibility stuff CHARACTER *3 MNTH(12) INTEGER *4 XLAST(4), YLAST(4), NLAST(4), NJUMP, NCOUNT INTEGER *4 N4HIST(32) INTEGER *4 N3HIST(32) REAL *4 NSUM(4), NNSUM(4), NUM8(4), NUM(4) REAL *8 VIS(4), VSUM(4), N2(4), VRMS(4) REAL *8 PHI, RPHI, PHIAVG, LASTPHI, PHIRMS, PHI12, LRPHI INCLUDE 'fritz.inc' C----------------------------------------------------------------------- REAL *8 JD CHARACTER *8 HMS EXTERNAL JD, HMS C----------------------------------------------------------------------- DATA MNTH / 'JAN', 'FEB', 'MAR', 'APR', 'MAY', 'JUN', $ 'JUL', 'AUG', 'SEP', 'OCT', 'NOV', 'DEC' / DATA FAST / .TRUE. / DATA N4HIST / 32*0 / DATA N3HIST / 32*0 / C----------------------------------------------------------------------- DRA = 0. DDEC = 0. DELAY2 = 0. DELAY3 = 0. C----------------------------------------------------------------------- C Print out headers, initialize unit numbers, read in the default C baselines. C CALL OPSTART ( versn, versd ) C----------------------------------------------------------------------- C Get and check the user-supplied input parameters; C open the input and output files; on error go back and C get new input parameters. MOREDATA is true if there C were no problems in starting up, STRTSTAR is the first C star to look at, STRTSEC and STOPSEC are the time limits C for the data in seconds from start of the year. C After all the data is read in, control returns here to C process another file, if desired. 1 CONTINUE CALL OPINPUT( MOREDATA, STRTYR, STRTMON, $ STRTDAY, STRTSTAR, STRTSEC, STOPSEC, $ INTTIM, MINAVG, ALL ) WRITE(TERMM,*) ' ' IF (.NOT. MOREDATA ) GO TO 900 TIME0 = DBLE (STRTSEC) ETIME = DBLE (STOPSEC) TDTOFF = 55 ISCAN = 0 WRITE(TERMP,*) ' Input file name = ', INDSN IF ( STRTMON .EQ. 0 ) THEN WRITE(TERMP,1020) ' Date = any ' ELSE WRITE(TERMP,1020) ' date = ', $ STRTYR, MNTH(STRTMON), STRTDAY END IF IF ( STRTSTAR .EQ. 0 ) THEN WRITE(TERMP,1030) ' Star = any ' ELSE WRITE(TERMP,1030) ' Star = ', STRTSTAR END IF WRITE(TERMP,*) ' Data averaged from ', HMS(STRTSEC), ' to ', $ HMS(STOPSEC) WRITE(TERMP,'(A,F6.3)') ' Integration time = ', INTTIM WRITE(TERMP,*) ' Minimum average = ', MINAVG, ' samples ' IF ( ALL ) THEN WRITE(TERMP,*) ' Write short averages to output file' ELSE WRITE(TERMP,*) ' Write scan averages to output file' END IF WRITE(TERMP,*) ' Data offset from start of file = ', OFFSET C----------------------------------------------------------------------- C Skip to start of the data to be processed. C 50 CONTINUE CALL RDATA ( NEWSTAR, MOREDATA, LOSTLOCK ) IF ( .NOT. MOREDATA ) GO TO 600 IF ( NEWSTAR ) GO TO 50 IF ((STRTYR .NE. IYEAR).AND.(STRTYR.NE.0)) THEN WRITE(TERMP,1050) 'year', STRTYR, 'year', IYEAR GO TO 600 ELSE IF ((STRTMON.NE.IMON).AND.(STRTMON.NE.0)) THEN WRITE(TERMP,1050) 'month', STRTMON, 'month', IMON GO TO 600 ELSE IF ((STRTDAY.NE.IDAY).AND.(STRTDAY.NE.0)) THEN WRITE(TERMP,1050) 'day', STRTDAY, 'day', IDAY GO TO 600 END IF IF ( MTIME .LT. TIME0 ) GO TO 50 IF ( MTIME .GT. ETIME ) GO TO 600 IF ( (ISTAR.NE.STRTSTAR).AND.(STRTSTAR.NE.0) ) GO TO 50 C TFIRST = MTIME TLAST = MTIME STRTYR = IYEAR STRTMON = IMON STRTDAY = IDAY WRITE(TERMP,1480) IDAY, MNTH(IMON), IYEAR C WRITE(TERMM,1490) IDAY, MNTH(IMON), IYEAR C DO 70 I = 1, NFILT C BDATA(I) = 0 C NLOCK = 0 C 70 CONTINUE C....................................................................... C We now have both time and the first data record. C One observation is a series of data on one star. DO 500 IOBS = 1, MAXOBS C C For each star, get C RAM, DECM the mean star position C PMRA, PMDEC, the proper motion C RV, PAR radial velocity and parallax C MAG, SPECT, the V magnitude, spectral class C ANGDIA the angular diameter in arcseconds CALL CATALOG( ISTAR, NAME, RAM, DECM, PMRA4, PMDEC4, $ PAR4, RV4, MAG4, BV, SPECT, ANGDIA ) PMRA = PMRA4 PMDEC = PMDEC4 PAR = PAR4 RV = RV4 MAG = MAG4 C....................................................................... C To speed up the averaging, the astrometry C need only be done once for each scan C C DJUL is the julian date at 0 hours UT in fractional days C UT1 is the universal time, in fractional days C TJD is the terrestial dynamical time, which differs from UT by C TDTOFF, the cummulative number of leap seconds (=55 on 1 Jan 1988) C C APSTAR computes the apparant place of the star (RAA, DECA) at the C current time corrected for everything except diurnal C abberation. C SIDTIM determines the sideral time in England (GAST) C DIURN computes the diurnal abberation and outputs the star C position to use for the rest of the calculation. C IF ( FAST ) THEN DJUL = JD ( IDAY, IMON, IYEAR ) UT1 = MTIME / 86400.D0 TJD = DJUL + UT1 + DBLE(TDTOFF)/86400.D0 NBODY = 3 CALL APSTAR (TJD, NBODY, RAM, DECM, PMRA, PMDEC, PAR, $ RV, RAA, DECA) CALL SIDTIM( DJUL, UT1, 1, GAST ) CALL DIURN ( TJD, GAST, RAA, DECA, RA, DEC ) END IF C C Twrite writes a time block to the visibility file C CALL TWRITE( INTTIM ) C....................................................................... C Choose the baseline for this observation BX = OBX (IBASE) BY = OBY (IBASE) BZ = OBZ (IBASE) CON= OCON(IBASE) C....................................................................... C Set up averages for each scan. C XSUM is the average time, relative to TO, the starting time. C YSUM is the average residual delay. C YYSYM is the average residual delay squared. C VSUM is the average visibility squared. C VRMS is the rms of the visibility squared. C NPTS is the number of averages in this scan. T0 = MTIME XSUM = 0.D0 YSUM = 0.D0 YYSUM = 0.D0 NPTS = 0.D0 DO 110 K = 1, NFILT VSUM(K) = 0.D0 VRMS(K) = 0.D0 110 CONTINUE DO 199 J=1,32 N4HIST(J) = 0 N3HIST(J) = 0 199 CONTINUE DO 300 I = 1, MAXCNT C----------------------------------------------------------------------- C Calculate the short averages, recommended times are C .5 seconds for one color data and visibilities C .1 seconds for two color data and central fringe id. C LASER0 = LASER TIME0 = MTIME ICNT8 = 0 DO J = 1, NFILT NUM (J) = 0. NSUM (J) = 0. NNSUM(J) = 0. NUM8 (J) = 0. END DO CALL LFIT ( 1, DELT, DLAS, ICNT, SLOPE,CONST,LRMS) C C Zero the accumulators needed for the phase calculations. C PHIRMS = 0. PHIAVG = 0. LRPHI = 0. LASTPHI = 0. PHI12 = 0. NJUMP = 0 NCOUNT = 0 DO 200 J = 1, 5000 DELT = MTIME- TIME0 DLAS = LASER- LASER0 CALL LFIT ( 2, DELT,DLAS,ICNT,SLOPE,CONST,LRMS) C C PHI is the observed phase in the tracking channel C RPHI is the unwrapped phase in the tracking channel. C C The unwrapping algorithm simply adds the appropriate multiple of C two pi to the observed phase to minimize the change in phase. C IF ( (X(1) .EQ. 0 ) .AND. ( Y(1) .EQ. 0 )) THEN PHI = 0. ELSE PHI = ATAN2 ( FLOAT(X(1)), FLOAT(Y(1)) ) END IF IF ( (LASTPHI-PHI) .GT. PI ) THEN NJUMP = NJUMP + 1 ELSE IF ( (LASTPHI-PHI) .LT. -PI ) THEN NJUMP = NJUMP - 1 END IF RPHI = PHI + 2.*PI*FLOAT(NJUMP) PHIAVG = PHIAVG + RPHI PHIRMS = PHIRMS + RPHI*RPHI PHI12 = PHI12 + LRPHI*RPHI LASTPHI = PHI LRPHI= RPHI NCOUNT = NCOUNT + 1 C C Accumulate the sums for the visibility average C DO K = 1, NFILT NUM (K) = NUM (K) + $ FLOAT( X(K)*X(K) + Y(K)*Y(K) - N(K) ) NSUM(K) = NSUM(K) + N(K) NNSUM(K) = NNSUM(K) + N(K)*N(K) END DO N4HIST(N(4)/16+1) = N4HIST(N(4)/16+1) + 1 N3HIST(N(3)/16+1) = N3HIST(N(3)/16+1) + 1 C C Do an 8 ms coherent average for visi data C IF ( MOD(J,2) .EQ. 1 ) THEN DO 160 K = 1, NFILT XLAST(K) = X(K) YLAST(K) = Y(K) NLAST(K) = N(K) 160 CONTINUE ELSE ICNT8 = ICNT8 + 1 DO 170 K = 1, NFILT NUM8(K) = NUM8(K) + FLOAT( $ ( X(K)+XLAST(K) ) **2 + $ ( Y(K)+YLAST(K) ) **2 - $ ( N(K)+NLAST(K) ) ) 170 CONTINUE END IF TLAST = MTIME CALL RDATA ( NEWSTAR, MOREDATA, LOSTLOCK ) C C Do no include data that spans a lost lock C IF ( .NOT. MOREDATA ) GO TO 201 IF ( NEWSTAR ) GO TO 201 IF ( LOSTLOCK ) GO TO 201 IF ( (MTIME-TIME0) .GE. (INTTIM-1.D-3)) GO TO 201 200 CONTINUE 201 CONTINUE IF ( ICNT .LT. MINAVG ) GO TO 299 C----------------------------------------------------------------------- C end of incoherent integration. CALL LFIT ( 3, DELT,DLAS,ICNT,SLOPE,CONST,LRMS) Comment WRITE(6,'(A,2F10.3)') ' LASER, RMS = ', DLAS, LRMS C C Modified 16 July 1990, David Mozurkewich. C Convert laser rms from microns to radians at tracking wavelength C (0.7 microns) so that it can be compared to the phase rms. C LRMS = LRMS * 2. * PI / 0.70 C TAVG = TIME0 + DELT LAVG = LASER0 + DLAS PHIAVG = PHIAVG / FLOAT(NCOUNT) PHIRMS = PHIRMS/FLOAT(NCOUNT) - PHIAVG*PHIAVG IF ( PHIRMS .LT. 0 ) THEN PHIRMS = - SQRT ( -PHIRMS ) ELSE PHIRMS = SQRT ( PHIRMS ) END IF C C Modified 3 July 90. Divided by NCOUNT-1 in next line, because C the first term in the sum is always zero. C PHI12 = PHI12 / FLOAT(NCOUNT-1) - PHIAVG*PHIAVG IF ( PHI12 .LT. 0. ) THEN PHI12 = - SQRT ( - PHI12 ) ELSE PHI12 = SQRT ( PHI12 ) END IF C DO 220 K = 1, NFILT NUM (K) = NUM (K) / FLOAT( ICNT ) NUM8 (K) = NUM8 (K) / FLOAT( ICNT8) NSUM (K) = NSUM (K) / FLOAT( ICNT ) NNSUM(K) = NNSUM(K) / FLOAT( ICNT ) IF ( NSUM(K) .EQ. 0. ) THEN Comment WRITE(TERMP,*) ' *** NSUM = 0 !!! channel ', K VIS(K) = 0. N2 (K) = 0. ELSE C C Here, VIS is the visibility (squared). VIS(K) = PI * PI * NUM (K) / $ ( 2.*NSUM(K)*NSUM(K) ) N2 (K) = N2(K) + NSUM(K) END IF VSUM(K) = VSUM(K) + VIS(K) VRMS(K) = VRMS(K) + VIS(K)*VIS(K) C C Now, VIS is the visibility COMMENT IF ( VIS(K) .LT. 0. ) THEN COMMENT VIS(K) = - SQRT( -VIS(K) ) COMMENT ELSE COMMENT VIS(K) = SQRT( VIS(K) ) COMMENT END IF 220 CONTINUE C....................................................................... C Determine the delay, based on the assumed baseline for the mean C time in this average. C DJUL = JD ( IDAY, IMON, IYEAR ) UT1 = TAVG / 86400.D0 TJD = DJUL + UT1 + DBLE(TDTOFF)/86400.D0 IF ( .NOT. FAST ) THEN NBODY = 3 CALL APSTAR (TJD, NBODY, RAM, DECM, PMRA, PMDEC, $ PAR, RV, RAA, DECA) CALL DIURN ( TJD, GAST, RAA, DECA, RA, DEC ) END IF CALL SIDTIM ( DJUL, UT1, 1, GAST ) HA = 15.D0 *( GAST - RA ) - LON IF ( HA .LT. -180.D0 ) HA = HA + 360.D0 DELAYC = CON + BX*COS(DEC/RAD)*COS(HA/RAD) $ - BY*COS(DEC/RAD)*SIN(HA/RAD) $ + BZ*SIN(DEC/RAD) C....................................................................... C Output to data files C Modified 3 July 1990. output LRMS instead of PHIRMS. C CALL VWRITE( TAVG, NUM, NUM8, NSUM, NNSUM, LRMS, PHI12) IF ( ALL ) THEN DELAY1 = LAVG WEIGHT = ICNT SIGMA = LRMS CALL WRITREC ( TERMO, UT1 ) END IF C C Average the difference delay. All of the data are C given the same weight. C NPTS = NPTS + 1.D0 IF ( I .EQ. 1 ) T0 = TAVG DTIME = DBLE( TAVG-T0 ) DLAS = LAVG - DELAYC YSUM = YSUM + DLAS XSUM = XSUM + DTIME YYSUM = YYSUM+ DLAS * DLAS 299 IF ( (.NOT. MOREDATA) .OR. NEWSTAR ) GO TO 301 300 CONTINUE 301 CONTINUE C IF ( NPTS .LT. 1.D0 ) THEN WRITE(TERMP,1570) IOBS, ISTAR, NAME, IBASE, $ MOD(IYEAR,100), IMON, IDAY GO TO 350 END IF XSUM = XSUM / NPTS YAVG = YSUM / NPTS LRMS = YYSUM / NPTS - YAVG * YAVG IF ( LRMS .GT. 0.D0 ) THEN LRMS = SQRT ( LRMS ) ELSE WRITE(6,*) ' laser rms is less than zero !!! ' LRMS = 0.D0 END IF WEIGHT= NPTS IF ( NPTS .EQ. 1 ) THEN DO 310 K = 1, NFILT VRMS(K) = 0. 310 CONTINUE ELSE DO 320 K = 1, NFILT N2 (K) = N2 (K) / NPTS VSUM(K) = VSUM(K) / NPTS VRMS(K) = VRMS(K)/NPTS - VSUM(K)*VSUM(K) IF ( VRMS(K) .LT. 0.D0 ) THEN VRMS(K) = 0. ELSE VRMS(K) = SQRT( VRMS(K) ) END IF COMMENT IF (VSUM(K) .LT. 0. ) THEN COMMENT VSUM(K) = - SQRT ( -VSUM(K) ) COMMENT IF ( ABS(VSUM(K)) .LT. .001 ) THEN COMMENT VRMS(K) = SQRT ( VRMS(K) ) COMMENT ELSE COMMENT VRMS(K) = - VRMS(K) / ( 2. * VSUM(K) ) COMMENT END IF COMMENT ELSE COMMENT VSUM(K) = SQRT ( VSUM(K) ) COMMENT IF ( ABS(VSUM(K)) .LT. .001 ) THEN COMMENT VRMS(K) = SQRT ( VRMS(K) ) COMMENT ELSE COMMENT VRMS(K) = VRMS(K) / ( 2. * VSUM(K) ) COMMENT END IF COMMENT END IF 320 CONTINUE END IF C....................................................................... C Calculated delay for the mean time DJUL = JD ( IDAY, IMON, IYEAR ) UT1 = ( T0 + XSUM ) / 86400.D0 CALL SIDTIM ( DJUL, UT1, 1, GAST ) TJD = DJUL + UT1 + DBLE(TDTOFF)/86400.D0 NBODY = 3 CALL APSTAR ( TJD, NBODY, RAM, DECM, PMRA, PMDEC, PAR, $ RV, RAA, DECA ) CALL DIURN ( TJD, GAST, RAA, DECA, RA, DEC ) HA = 15.D0 *( GAST - RA ) - LON IF ( HA .LT. -180.D0 ) HA = HA + 360.D0 C....................................................................... DELAYC = CON + BX*COS(DEC/RAD)*COS(HA/RAD) $ - BY*COS(DEC/RAD)*SIN(HA/RAD) $ + BZ*SIN(DEC/RAD) ZA = RAD * ACOS( SIN(DEC/RAD)*SIN(LAT/RAD) + $ COS(DEC/RAD)*COS(LAT/RAD)*COS(HA/RAD) ) TIME0 = DJUL DELAY1 = DELAYC + YAVG DELAY2 = DELAYC DELAY3 = YAVG C....................................................................... SIGMA = LRMS IF ( .NOT. ALL ) THEN CALL WRITREC ( TERMO, UT1 ) END IF TEMP = INT ( 8.64D4*UT1 ) DO 345 K = 1, 4 FDARK(K) = DARK(K)/1250. 345 CONTINUE IF (MAKEHI.GT.0) WRITE(TERMP,1551) N3HIST IF (MAKEHI.GT.0) WRITE(TERMP,1551) N4HIST 1551 FORMAT (32I5) WRITE(TERMP,1550) IOBS, ISTAR, NAME, IBASE, $ MOD(IYEAR,100), IMON, IDAY, HMS(TEMP), $ NLOCK, (4*NDATA)/NLOCK, DELAY3, $ LRMS, (VSUM(K), K = 1, 4), BDATA WRITE(TERMM,1590) IOBS 350 CONTINUE C C Write dark data to visibilities files. Note that the dark counts C are set = 0 at the start of each scan, so if there is no dark C data for this scan, none is output. The dark count record is needed C in the next stage of processing for the other stuff in the record. C HA0 = HA DEC0 = DEC CALL DWRITE( UT1, DELAY1, TFIRST, TLAST, NLOCK ) C IF ( .NOT. MOREDATA ) GO TO 600 IF ( MTIME .GT. ETIME ) GO TO 600 IF ( NEWSTAR ) THEN NLOCK = 1 NDATA = 0 DO K = 1, NFILT BDATA(K) = 0 NLOCK = 0 END DO CALL RDATA ( NEWSTAR, MOREDATA, LOSTLOCK ) IF ( .NOT. MOREDATA ) GO TO 600 IF ( NEWSTAR ) THEN WRITE(TERMP,*) ' Multiple new stars ' GO TO 350 END IF IF ( .NOT. MOREDATA ) GO TO 600 IF ( MTIME .GT. ETIME ) GO TO 600 TFIRST = MTIME TLAST = MTIME END IF C C If strtstar is set, check to see if this is the correct star. C IF ( STRTSTAR .NE. 0 ) THEN 450 CONTINUE IF ( STRTSTAR .NE. ISTAR ) THEN CALL RDATA ( NEWSTAR, MOREDATA, LOSTLOCK ) IF ( .NOT. MOREDATA ) GO TO 600 IF ( MTIME .GT. ETIME ) GO TO 600 IF ( NEWSTAR ) GO TO 350 GO TO 450 END IF END IF 500 CONTINUE 600 CONTINUE C Finish up and exit. WRITE(TERMP,1610) IREC CALL OPCLOSE GO TO 1 900 CONTINUE STOP 1010 FORMAT ( A, 2F11.0 ) 1020 FORMAT ( A, I5, A4, I5 ) 1030 FORMAT ( A, I5 ) 1050 FORMAT(' starting ', A, I5, ' is not equal to current ', A, I5) 1400 FORMAT ( I5, I2, F15.3, 4F8.1, 4F7.3 ) 1480 FORMAT ( ' Averaged delays for ', I5, A4, I5, 9X, / $ ' star numb m.sec', $ ' resid rms VISIBILITY SQUARED bad data' / $ ' No. star name B YY MM DD HH:MM:SS locks lock', $ ' laser las CH1 CH2 CH3 CH4 b1 b2 b3 b4' ) 1490 FORMAT ( ' DATE = ', I5, A4, I5 / $ ' star numb m.sec', $ ' bad data resid' / $ ' No. star B name HH:MM:SS locks lock', $ ' b1 b2 b3 b4 laser' ) 1550 FORMAT ( I4, I5, 1X, '"', A8, '"', I3, 3(1X, I2.2), 1X, A8, $ I6, I7, F8.1, $ F10.1, 4F5.2, 1X, 4I3 ) 1570 FORMAT ( I4, I5, 1X, '"', A8, '"', I3, 3(1X, I2.2), 5X, $ ' NO DATA ' ) 1590 FORMAT ( '+ FINISHED SCAN NUMBER ', I5, $ ) 1610 FORMAT ( ' FRITZ ends. Block number ', I5 ) END SUBROUTINE LFIT( MODE, X0, Y0, N, SLOPE, CONST, SIGMA ) C C Least squares fit of a line through the points, x, y. IMPLICIT NONE INTEGER *4 N, MODE REAL *8 X0, Y0, SLOPE, CONST, SIGMA REAL *8 X, Y, XX, YY, XY IF ( MODE .EQ. 1 ) THEN C INITIALIZE X = 0. Y = 0. XY = 0. XX = 0. YY = 0. N = 0 ELSE IF ( MODE .EQ. 2 ) THEN C INCREMENT X = X + X0 XX = XX + X0 * X0 XY = XY + X0 * Y0 YY = YY + Y0 * Y0 Y = Y + Y0 N = N + 1 ELSE IF ( (MODE.EQ.3).AND.(N.LT.5) ) THEN Comment WRITE(0,*) ' NOT ENOUGH POINTS TO SOLVE ' SIGMA = 0.D0 ELSE IF ( MODE .EQ. 3 ) THEN C SOLVE X0 = X / FLOAT(N) Y0 = Y / FLOAT(N) XY = XY / FLOAT(N) XX = XX / FLOAT(N) YY = YY / FLOAT(N) SLOPE = ( XY - X0*Y0 ) / ( XX - X0*X0 ) CONST = Y0 - SLOPE*X0 COMMENT SIGMA = YY - Y0*Y0 + SLOPE*( XY -X0*X0 - 2.*XY + 2.*X0*Y0 ) SIGMA = YY + SLOPE*SLOPE*XX + CONST**2 - 2.*SLOPE*XY $ - 2.*CONST*Y0 + 2.*CONST*SLOPE*X0 SIGMA = SQRT ( MAX(0., SIGMA) ) ELSE WRITE(0,*) ' VALUE OF MODE ', MODE, ' IS NOT ALLOWED ' END IF RETURN END