SUBROUTINE DARKPLOT C C This subroutine plots dark counts versus time or hour angle C fits a polynomial to the data, allows the user to delete C points and and will make hard copies on request. C IMPLICIT UNDEFINED(A-Z) SAVE INCLUDE 'VPLOT.INC' INTEGER *4 I, J, SYMBOL(5), ICHAN, IPLOT, ICUR, IFRAME INTEGER *4 POINTER(CNT), MASK(4), OKDARK(CNT), LIST(CNT) INTEGER *4 IERR, NDF REAL *4 XMIN, XMAX, YMIN, YMAX, RX, RY, DX, DY REAL *4 XLO(4), XHI(4), YLO(4), YHI(4), BARX(2), BARY(2), DELTA REAL *4 DIST, DIST2, XPLOT(CNT,4), YPLOT(CNT,4), YFIT(CNT,4) REAL *4 ERROR(CNT), SCALE, XLINE(CNT), YLINE(CNT) REAL *4 WIDTH, XPORT(2), YPORT(2), XMID, YMID REAL *8 COEF(CNT) CHARACTER *60 TITLE, INS CHARACTER * 1 CH LOGICAL HCOPY INCLUDE 'HARDCOPY.INC' C----------------------------------------------------------------------- DATA SYMBOL / 16, 3, 5, 6, 13 / DATA ICHAN, ICUR, IFRAME, D_ORDER / 1, 1, 1, 8*0 / DATA MASK / 1, 2, 4, 8 / DATA XPORT, YPORT / .20, .80, .20, .80 / HCOPY = .FALSE. C----------------------------------------------------------------------- IERR = 1 DO WHILE ( ((CH.NE.'N').AND.(CH.NE.'Y')) .OR. (IERR.NE.0) ) WRITE(6,*) ' Restore deleted scans (Y/N)?' READ (5,'(A)',IOSTAT=IERR) CH CALL CAPS ( CH ) END DO IF ( CH .EQ. 'Y' ) THEN DO I = 1, ISCAN C C Do not undelete a scan if there is no data. C IF ( HOURS(I) .EQ. 0.0 ) THEN OKDARK(I) = 15 ELSE OKDARK(I) = 0 END IF END DO ELSE DO I = 1, ISCAN OKDARK(I) = OKSCAN(I) END DO END IF C----------------------------------------------------------------------- C Dusk and Dawn are defined as the first and last points included C in the dark count fits. The data must be sorted by time or else C Dawn and Dusk will not be calculated correctly. C IF ( (IDUSK .EQ. 0).OR.(IDAWN .EQ. 0) ) THEN IDUSK = 1 IDAWN = ISCAN ELSE 190 CONTINUE WRITE(6,'(A,$)') ' use old dusk and dawn values (Y/N) ?' READ (5,*) CH CALL CAPS (CH) IF ( CH .EQ. 'Y' ) THEN DO I = 1, IDUSK-1 OKDARK(I) = 15 END DO DO I = 1, IDUSK-1 OKDARK(I) = 15 END DO ELSE IF ( CH .NE. 'N' ) THEN GO TO 190 END IF END IF C---------------------------------------------------------------------- C Return here for start/end of hard copy, otherwise loop starts at 210. C 200 CONTINUE IF ( HCOPY ) THEN NPLOTS= NPLOTS + 1 WRITE(PLOTFILE(6:8),'(I3.3)') NPLOTS CALL PGBEGIN(0, PLOTFILE//PLOT_DEV, 1, 1 ) WRITE(6,*) ' CREATING FILE ', PLOTFILE ELSE CALL PGBEGIN(0, '/EGA', 1, 1) END IF 210 CONTINUE CALL PGADVANCE IPLOT = 0 C C Get data for this channel C XMIN = 1.E36 XMAX = -1.E36 YMIN = 1.E36 YMAX = -1.E36 DO I = 1, ISCAN IF ( ( D(I,ICHAN) .GT. 0. ) .AND. $ ( IAND(OKDARK(I),MASK(ICHAN)) .EQ. 0 ) ) THEN IPLOT = IPLOT + 1 XPLOT(IPLOT,1) = HOURS(I) XPLOT(IPLOT,2) = DELAY(I) XPLOT(IPLOT,3) = HOURS(I) XPLOT(IPLOT,4) = DELAY(I) YPLOT(IPLOT,1) = D(I,ICHAN) YPLOT(IPLOT,2) = D(I,ICHAN) C C Dark counts are measured for 5 seconds. The value in D is for one C coherent integration. The next line fixed, so that it works for C any coherent integration time. (4 Dec 1990 DM) C ERROR(IPLOT) = SQRT ( D(I,ICHAN) * FLOAT(COHINT) / 5000. ) POINTER(IPLOT) = I C XMIN = MIN ( XPLOT(1,1), XMIN ) XMAX = MAX ( XPLOT(1,1), XMAX ) YMIN = MIN ( YPLOT(1,1), YMIN ) YMAX = MAX ( YPLOT(1,1), YMAX ) XMIN = MIN ( XPLOT(1,2), XMIN ) XMAX = MAX ( XPLOT(1,2), XMAX ) END IF END DO IF ( IPLOT .EQ. 0 ) THEN D_ORDER(1,ICHAN) = 0 D_ORDER(2,ICHAN) = 0 D_RMS(ICHAN) = 0. D_CHISQ(ICHAN) = 0. GO TO 300 END IF IDUSK = MAX(IDUSK, POINTER(1)) IDAWN = MIN(IDAWN, POINTER(IPLOT)) C----------------------------------------------------------------------- C Do a fit to background versus time. C IF ( D_ORDER(1,ICHAN) .LT.0 ) THEN WIDTH = - D_ORDER(1,ICHAN) CALL BOXCAR ( XPLOT(1,1), YPLOT(1,1), YFIT(1,1), IPLOT, $ XMIN, XMAX, YMIN, YMAX, WIDTH, $ HOURS, DARKFIT(1,ICHAN), ISCAN, D_RMS(ICHAN)) ELSE CALL FIT_1D ( XPLOT(1,1), YPLOT(1,1), YFIT(1,1), IPLOT, $ D_ORDER(1,ICHAN), D_RMS(ICHAN), COEF ) END IF C----------------------------------------------------------------------- C Do a fit to background versus delay. C CALL FIT_1D ( XPLOT(1,2), YPLOT(1,2), YFIT(1,2), IPLOT, $ D_ORDER(2,ICHAN), D_RMS(ICHAN), COEF ) C----------------------------------------------------------------------- C Do the combined fit. The fit is simultaneous unless the C fit to time is boxcar smoothing ( ORDER < 0 ) C DARKFIT contains the fit, for all scans. It is used for editing. C IF ( D_ORDER(1,ICHAN) .LT. 0 ) THEN DO I = 1, IPLOT YPLOT(I,3) = YPLOT(I,1) - YFIT(I,1) END DO CALL FIT_1D ( XPLOT(1,2), YPLOT(1,3), YFIT(1,3), IPLOT, $ D_ORDER(2,ICHAN), D_RMS(ICHAN), COEF ) DO I = 1, IPLOT YPLOT(I,3) = YPLOT(I,3) - YFIT(I,3) YPLOT(I,4) = YPLOT(I,3) END DO DO I = 1, ISCAN DARKFIT(I,ICHAN) = DARKFIT(I,ICHAN) + COEF(1) DO J = 1, D_ORDER(2,ICHAN) DARKFIT(I,ICHAN) = DARKFIT(I,ICHAN) $ + COEF(J+1) * ( DELAY(I)**J ) END DO END DO ELSE CALL FIT_2D ( XPLOT(1,1), YPLOT(1,1), YFIT(1,3), IPLOT, $ D_ORDER(1,ICHAN), D_RMS(ICHAN), CNT, COEF ) DO I = 1, IPLOT Comment DARKFIT(I,ICHAN) = YFIT(I,3) YPLOT(I,3) = YPLOT(I,1) - YFIT(I,3) YPLOT(I,4) = YPLOT(I,3) END DO DO I = 1, ISCAN DARKFIT(I,ICHAN) = COEF(1) DO J = 1, D_ORDER(1,ICHAN) DARKFIT(I,ICHAN) = $ DARKFIT(I,ICHAN) + COEF(J+1) * ( HOURS(I)**J ) END DO DO J = 1, D_ORDER(2,ICHAN) DARKFIT(I,ICHAN) = DARKFIT(I,ICHAN) $ + COEF(J+1+D_ORDER(1,ICHAN)) * ( DELAY(I)**J ) END DO END DO END IF C----------------------------------------------------------------------- C Evaluate the range of fitted values. C D_RANGE(1,ICHAN) = 100.0 D_RANGE(2,ICHAN) = 0.0 DO I = 1, ISCAN D_RANGE(1,ICHAN) = MIN ( D_RANGE(1,ICHAN), DARKFIT(I,ICHAN) ) D_RANGE(2,ICHAN) = MAX ( D_RANGE(2,ICHAN), DARKFIT(I,ICHAN) ) END DO C----------------------------------------------------------------------- C Evaluate the reduced chi square. C D_CHISQ(ICHAN) = 0. IF ( D_ORDER(1,ICHAN) .LT. 0 ) THEN NDF = -60.*(HOURS(IDAWN)-HOURS(IDUSK))/D_ORDER(1,ICHAN) ELSE NDF = D_ORDER(1,ICHAN) END IF NDF = IPLOT - NDF - D_ORDER(2,ICHAN) - 1 IF ( NDF .LE. 0 ) THEN D_CHISQ(ICHAN) = 999.99 ELSE DO I = 1, IPLOT DELTA = YPLOT(I,3) / ERROR(I) D_CHISQ(ICHAN) = D_CHISQ(ICHAN) + DELTA*DELTA END DO IF ( IPLOT .GT. 1 ) THEN D_CHISQ(ICHAN) = D_CHISQ(ICHAN) / FLOAT(NDF) END IF END IF C----------------------------------------------------------------------- C Sort the data by delay so that the plots will look nice. CALL SORT ( IPLOT, XPLOT(1,2), LIST ) DO I = 1, IPLOT XLINE(I) = XPLOT(LIST(I),2) YLINE(I) = YFIT (LIST(I),2) END DO C----------------------------------------------------------------------- C Make 4 plots: dark and redisual dark versus time and delay. C 300 CONTINUE YMID = 0.5 * ( YPORT(1) + YPORT(2) ) XMID = 0.5 * ( XPORT(1) + XPORT(2) ) CALL PGSCI ( 5 ) CALL PGSCH ( 1.5 ) CALL PGVPORT ( XPORT(1), XPORT(2), YPORT(1), YPORT(2) ) WRITE(TITLE, 1305) ICHAN, FILTERID(ICHAN) CALL PGSCH ( 1.0 ) CALL PGMTEXT ( 'T', 2.0, 0.5, 0.5, TITLE ) WRITE(TITLE, 1306) D_ORDER(1,ICHAN), D_ORDER(2,ICHAN) CALL PGMTEXT ( 'R', 5.0, 0.5, 0.5, TITLE ) WRITE(TITLE, 1310) D_RMS(ICHAN), D_CHISQ(ICHAN) CALL PGMTEXT ( 'B', 5.0, 0.5, 0.5, TITLE ) C DO J = 1, 4 IF ( IPLOT .GT. 0 ) THEN XMIN = XPLOT(1,J) XMAX = XPLOT(1,J) YMIN = YPLOT(1,J) YMAX = YPLOT(1,J) DO I = 2, IPLOT XMIN = MIN ( XPLOT(I,J), XMIN ) XMAX = MAX ( XPLOT(I,J), XMAX ) YMIN = MIN ( YPLOT(I,J), YMIN ) YMAX = MAX ( YPLOT(I,J), YMAX ) END DO ELSE XMIN = 0. XMAX = 10. YMIN = 0. YMAX = 1. END IF C Set big letters so that PGENV will leave C more room around the plot. CALL PGSCI ( 5 ) CALL PGSLW ( 1 ) CALL PGSCH ( 1.) CALL PGRNGE ( XMIN, XMAX, XLO(J), XHI(J) ) CALL PGRNGE ( YMIN, YMAX, YLO(J), YHI(J) ) IF ( J .EQ. 1 ) THEN CALL PGVPORT ( XPORT(1), XMID, YMID, YPORT(2) ) TITLE = 'DARK COUNT / 4 MS' CALL PGMTEXT ( 'L', 3.5, 0.5, 0.5, TITLE ) CALL PGWINDOW ( XLO(J), XHI(J), YLO(J), YHI(J) ) CALL PGBOX ( 'BCMST', 0.0, 0, 'BCNVST', 0.0, 0 ) ELSE IF ( J .EQ. 2 ) THEN CALL PGVPORT ( XMID, XPORT(2), YMID, YPORT(2) ) CALL PGWINDOW ( XLO(J), XHI(J), YLO(J), YHI(J) ) CALL PGBOX ( 'BCMST', 0.0, 0, 'BCMVST', 0.0, 0 ) ELSE IF ( J .EQ. 3 ) THEN CALL PGVPORT ( XPORT(1), XMID, YPORT(1), YMID ) TITLE = 'RESIDUAL DARK COUNT' CALL PGMTEXT ( 'L', 3.5, 0.5, 0.5, TITLE ) TITLE = 'TIME' CALL PGMTEXT ( 'B', 1.8, 0.5, 0.5, TITLE ) CALL PGWINDOW ( XLO(J), XHI(J), YLO(J), YHI(J) ) CALL PGBOX ( 'ABCNST', 0.0, 0, 'BCNVST', 0.0, 0 ) ELSE IF ( J .EQ. 4 ) THEN CALL PGVPORT ( XMID, XPORT(2), YPORT(1), YMID ) TITLE = 'DELAY' CALL PGMTEXT ( 'B', 1.8, 0.5, 0.5, TITLE ) CALL PGWINDOW ( XLO(J), XHI(J), YLO(J), YHI(J) ) CALL PGBOX ( 'ABCNST', 0.0, 0, 'BCMVST', 0.0, 0 ) END IF C CALL PGSCI ( 6 ) DO I = 1, IPLOT BARX(1) = XPLOT(I,J) BARX(2) = XPLOT(I,J) BARY(1) = YPLOT(I,J) - ERROR(I) BARY(2) = YPLOT(I,J) + ERROR(I) CALL PGLINE( 2, BARX, BARY ) END DO C----------------------------------------------------------------------- C Plot the fitted line through the data. C IF ( IPLOT .GT. 2 ) THEN IF ( J .EQ. 1 ) THEN CALL PGLINE ( IPLOT, XPLOT(1,J), YFIT(1,J) ) ELSE IF ( J .EQ. 2 ) THEN CALL PGLINE ( IPLOT, XLINE, YLINE ) END IF END IF END DO IF ( HCOPY ) THEN CALL PGEND HCOPY = .FALSE. GO TO 200 END IF C C The plot is plotted, now ask for command C Switch to text mode C CALL PGVPORT ( XPORT(1), XPORT(2), YPORT(1), YPORT(2) ) CALL PGWINDOW ( 0.0, 2.0, 0.0, 2.0 ) CALL EGA_RESTORE_DEFAULT 320 CONTINUE RX = ( XPLOT(ICUR,IFRAME) - XLO(IFRAME) ) $ / ( XHI(IFRAME) - XLO(IFRAME) ) RY = ( YPLOT(ICUR,IFRAME) - YLO(IFRAME) ) $ / ( YHI(IFRAME) - YLO(IFRAME) ) IF ( (IFRAME .EQ. 1) .OR. (IFRAME .EQ. 2) ) THEN RY = RY + 1.0 END IF IF ( (IFRAME .EQ. 2) .OR. (IFRAME .EQ. 4) ) THEN RX = RX + 1.0 END IF CALL LOCATE(0,0) CALL WRITE_STRING ( 'Command:' ) CALL PGCURSE(RX, RY, CH ) CALL CAPS ( CH ) IF ( CH .EQ. '?' ) THEN CALL LOCATE (0, 1) CALL WRITE_STRING ( ' ? Displays this ' ) CALL LOCATE (0, 2) CALL WRITE_STRING ( ' C toggle channel ' ) CALL LOCATE (0, 3) CALL WRITE_STRING ( ' H make hardcopy ' ) CALL LOCATE (0, 4) CALL WRITE_STRING ( ' I identify a scan ' ) CALL LOCATE (0, 5) CALL WRITE_STRING ( ' D delete a point ' ) CALL LOCATE (0, 6) CALL WRITE_STRING ( ' O change order of fit ' ) CALL LOCATE (0, 7) CALL WRITE_STRING ( ' R move right one point' ) CALL LOCATE (0, 8) CALL WRITE_STRING ( ' L move left one point' ) CALL LOCATE (0, 9) CALL WRITE_STRING ( ' T change frames ' ) CALL LOCATE (0,10) CALL WRITE_STRING ( ' Q quit ' ) CALL LOCATE (0,11) GO TO 320 ELSE IF ( CH .EQ. 'C' ) THEN C Toggle channel ICHAN = MOD ( ICHAN, 4 ) + 1 ELSE IF ( CH .EQ. 'H' ) THEN C Make a hard copy HCOPY = .TRUE. CALL PGEND GO TO 200 ELSE IF ( CH .EQ. 'I' ) THEN C C Identify the cursor point C CALL LOCATE ( 0, 1 ) WRITE(6,*) ' STAR = ', STAR(POINTER(ICUR)), $ ' SCAN = ', POINTER(ICUR) GO TO 320 ELSE IF ( CH .EQ. 'D' ) THEN C C Delete the current point (but only from the dark count fitting). C OKDARK(POINTER(ICUR)) = $ IOR( OKDARK(POINTER(ICUR)),MASK(ICHAN)) IF ( ICUR .EQ. IPLOT ) ICUR = IPLOT - 1 C change the order of the fit ELSE IF ( CH .EQ. 'O' ) THEN CALL LOCATE (0, 0 ) WRITE(INS,1400) D_ORDER(1,ICHAN), D_ORDER(2,ICHAN) CALL WRITE_STRING( INS ) READ( 5,*) D_ORDER(1,ICHAN), D_ORDER(2,ICHAN) C C Check limits C IF ( D_ORDER(1,ICHAN) .GE. 0 ) THEN D_ORDER(1,ICHAN) = MIN( 10, D_ORDER(1,ICHAN) ) ELSE D_ORDER(1,ICHAN) = MIN ( -25., D_ORDER(1,ICHAN) ) END IF D_ORDER(2,ICHAN) = MIN( 10, MAX( D_ORDER(2,ICHAN), 0 ) ) C ELSE IF ( CH .EQ. 'R' ) THEN ICUR = MOD( ICUR, IPLOT ) + 1 GO TO 320 C ELSE IF ( CH .EQ. 'L' ) THEN ICUR = MOD( ICUR-2+IPLOT, IPLOT ) + 1 GO TO 320 C ELSE IF ( CH .EQ. 'T' ) THEN IFRAME = MOD( IFRAME, 4 ) + 1 GO TO 320 C quit ELSE IF ( CH .EQ. 'Q' ) THEN GO TO 900 ELSE WRITE(6,*) ' INVALID COMMAND ' END IF CALL PGASK ( .FALSE. ) CALL PGADVANCE GO TO 210 900 CONTINUE CALL PGEND WRITE(6,*) ' Do you want to use these dark count curves ' WRITE(6,*) ' to edit the dark count data (Y/N)?' READ (5,'(A)') CH CALL CAPS ( CH ) IF ( CH .EQ. 'Y' ) THEN CALL DARKEDIT ( OKDARK ) ELSE IF ( CH .NE. 'N' ) THEN GO TO 900 END IF RETURN 1305 FORMAT ( 'CHANNEL', I3, ' = ', A7, 20X ) 1306 FORMAT ( 'ORDER OF FITS TIME =', I5, ' DELAY =', I4 ) 1310 FORMAT ( 'RMS OF FIT =', F8.3, ' CHISQ/DOF =', F8.1 ) 1400 FORMAT ( ' ORDER OF FITS ARE', 2I4, ' INPUT NEW VALUES:' ) END SUBROUTINE DARKEDIT ( OKDARK ) C C Does automatic editing of the dark count data. C C VERSION 2.0 13 April 1989 C C This program uses the estimated dark counts and dark count errors, C generated by DARKPLOT. C C The user is requested to input a minimum allowed uncertainty. If C points outside of twilight deviate by more than this fraction of C of the star's intensity, and this is more than NSIG*the expected C rms for this star, then the data is deleted. C C A table of rms of the dark count, normalized to the strength of the C star is displayed. C IMPLICIT UNDEFINED (A-Z) SAVE INCLUDE 'VPLOT.INC' INCLUDE 'HARDCOPY.INC' INTEGER *4 SYMBOL(5), SPTR, N, I, J, K, L, NPTS(4) INTEGER *4 ISTAR, ICHAN, MASK(4), WHERE(CNT), OKDARK(*) REAL *4 XPLOT(CNT), YPLOT(CNT,4), ZPLOT(CNT,4) REAL *4 YFIT(CNT), RMS(4), DELT, NSIG, F REAL *4 RX, RY, DAWN, DUSK, DELTA REAL *4 DAVG(4), DRMS(4), FLUX(4), XPORT(2), YPORT(2) REAL *4 XMIN, XMAX, YMIN(4), YMAX(4), XHI, XLO, YHI, YLO LOGICAL *4 HCOPY, FIXFLAG CHARACTER * 1 CH CHARACTER *40 TITLE C********************************************************************** DATA SYMBOL / 16, 3, 5, 6, 13 / DATA MASK / 1, 2, 4, 8 / DATA XPORT, YPORT / .20, .80, .10, .90 / C---------------------------------------------------------------------- C Times of dawn and dusk are determined in darkplot. They are no C longer determined here. DUSK = HOURS(IDUSK) DAWN = HOURS(IDAWN) C---------------------------------------------------------------------- DO ICHAN = 1, 4 RMS(ICHAN) = 0. END DO N = 0 DO 90 I = IDUSK, IDAWN DO ICHAN = 1, 4 IF ( D(I,ICHAN) .EQ. 0.0 ) GO TO 90 END DO N = N + 1 DO ICHAN = 1, 4 DELTA = D(I,ICHAN) - DARKFIT(I,ICHAN) RMS(ICHAN) = RMS(ICHAN) + DELTA * DELTA END DO 90 CONTINUE DO ICHAN = 1, 4 RMS(ICHAN) = SQRT ( RMS(ICHAN) / FLOAT(N) ) END DO C----------------------------------------------------------------------- WRITE(6,1100) ' DUSK = ', DUSK, ' hours point ', IDUSK WRITE(6,1100) ' DAWN = ', DAWN, ' hours point ', IDAWN WRITE(6,1110) ' RMS OF BACKGROUND = ', RMS, ' counts/4ms ' WRITE(6,*) ' If the measured background for a scan deviates ' WRITE(6,*) ' from the actual value by more than a fraction ' WRITE(6,*) ' of the number of photons from the star, then ' WRITE(6,*) ' the observed visibility will also be wrong by ' WRITE(6,*) ' by about the same fraction. Scans with deviant ' WRITE(6,*) ' backgrounds will be flagged ' WRITE(6,*) ' Input maximum allowed deviation ' 100 CONTINUE WRITE(6,'(A,$)') ' Recommended value = .03:' READ (5,*) F IF ( ( F .GT. .5 ) .OR. ( F .LT. .0005 ) ) THEN WRITE(6,*) ' YOU''VE GOT TO BE JOKING' GO TO 100 END IF 110 CONTINUE WRITE(6,*) ' Descrepant twilight backgrounds are always deleted' WRITE(6,*) ' Delete or fix other descrepant scans (D/F)?' READ (5,'(A)') CH CALL CAPS ( CH ) IF ( CH .EQ. 'F' ) THEN FIXFLAG = .TRUE. ELSE IF ( CH .EQ. 'D' ) THEN FIXFLAG = .FALSE. ELSE GO TO 110 END IF C----------------------------------------------------------------------- HCOPY = .FALSE. ICHAN = 1 NSIG = 3.0 CALL PGBEGIN(0, '/EGA', 1, 1) CALL PGADVANCE 1100 FORMAT ( A, F10.2, A ) 1110 FORMAT ( A, 4F10.2 ) 1500 FORMAT ( 'CH =', I2, A ) 1510 FORMAT ( ' STAR = ', I5, 2X, A8 ) 1530 FORMAT ( ' = ', F7.5, ' FIT RMS = ', F7.5 ) C---------------------------------------------------------------------- C Plot the dark counts for each star. C C For each star/channel pair, determine the mean offset between C the measured background and the fitted background. Stars close to C the moon, or stars in the brighter, southern half of the sky often C have a mean background somewhat higher than the fitted value. The C offset is determined using the scans that were not deleted in the C dark count fitting. Otherwise one scan with a spurious dark count C can corrupt the offset determined for that star. If an entire star C was deleted from the dark count fitting, the offset is set equal to C zero for that star. C C Missing data for a star should be given a value that agrees with C the star data in amplitude, but agrees with the general trends C in the fit. C C DARKAVG is the mean background for that star. C DARKRMS is the rms of the measured background for that star C around the fitted value. C DARKOFF is the mean difference between the measured background C and the fitted background. C C FLUX is the photon rate observed in each channel. The scan should C not be used if a large fraction of the observed photons come from C background. C XMIN = 23.0 XMAX = 0.0 DO I = 1, ISCAN IF ( HOURS(I) .NE. 0.D0 ) THEN XMIN = MIN ( XMIN, HOURS(I) ) XMAX = MAX ( XMAX, HOURS(I) ) END IF END DO DO 350 SPTR = 1, NSTAR N = 0 DO ICHAN = 1, 4 YMIN(ICHAN) = 0. YMAX(ICHAN) = 0. FLUX(ICHAN) = 0. DAVG(ICHAN) = 0. DRMS(ICHAN) = 0. DARKAVG(SPTR,ICHAN) = 0. DARKRMS(SPTR,ICHAN) = 0. DARKOFF(SPTR,ICHAN) = 0. NPTS(ICHAN) = 0 END DO DO 220 I = 1, ISCAN IF ( HOURS(I) .EQ. 0.D0 ) GO TO 220 IF ( STAR(I) .NE. SLIST(SPTR) ) GO TO 220 N = N + 1 WHERE(N) = I XPLOT(N) = HOURS(I) DO ICHAN = 1, 4 YPLOT(N,ICHAN) = D(I,ICHAN) ZPLOT(N,ICHAN) = DARKFIT(I,ICHAN) YMAX(ICHAN) = MAX( YMAX(ICHAN), YPLOT(N,ICHAN) ) C C If there is dark data and flux data for this scan and if it was C used in the dark data fit, then it should be used here. C IF ( ( D(I,ICHAN) .NE. 0. ) .AND. $ ( IAND(OKSCAN(I),MASK(ICHAN)) .EQ. 0 ) .AND. $ ( IAND(OKDARK(I),MASK(ICHAN)) .EQ. 0 ) ) THEN FLUX(ICHAN) = FLUX(ICHAN) + BDEN(I,ICHAN) DELTA = YPLOT(N,ICHAN) - ZPLOT(N,ICHAN) DARKAVG(SPTR,ICHAN) = DARKAVG(SPTR,ICHAN) $ + YPLOT(N,ICHAN) DARKRMS(SPTR,ICHAN) = DARKRMS(SPTR,ICHAN) $ + DELTA * DELTA DARKOFF(SPTR,ICHAN) = DARKOFF(SPTR,ICHAN) $ + DELTA NPTS(ICHAN) = NPTS(ICHAN) + 1 DAVG(ICHAN) = DAVG(ICHAN) + YPLOT(N,ICHAN) DRMS(ICHAN) = DRMS(ICHAN) $ + YPLOT(N,ICHAN) * YPLOT(N,ICHAN) END IF END DO 220 CONTINUE C---------------------------------------------------------------------- C Make four plots per page, one for each channel. C CALL PGSCI ( 5 ) CALL PGVPORT ( XPORT(1), XPORT(2), YPORT(1), YPORT(2) ) CALL PGSLS ( 1 ) CALL PGSCH ( 1.5 ) TITLE = 'TIME (HOURS)' CALL PGMTEXT ( 'B', 1.8, 0.5, 0.5, TITLE ) TITLE = 'BACKGROUND COUNT' CALL PGMTEXT ( 'L', 3.5, 0.5, 0.5, TITLE ) WRITE(TITLE,1510) SLIST(SPTR), STARNAME(SPTR) CALL PGMTEXT ( 'T', 2.0, 0.5, 0.5, TITLE ) TITLE = 'BACKGROUND/STAR FLUX' CALL PGMTEXT ( 'R', 4.5, 0.5, 0.5, TITLE ) CALL PGSCH ( 1.0 ) TITLE = ' ' CALL PGMTEXT ( 'T', 1.9, 0.5, 0.5, TITLE ) TITLE = ' ' CALL PGMTEXT ( 'B', 4.0, 0.5, 0.5, TITLE ) CALL PGSLW ( 1 ) DO ICHAN = 1, 4 IF ( NPTS(ICHAN) .EQ. 0 ) THEN DRMS(ICHAN) = 0. DARKRMS(SPTR,ICHAN) = 0. YMAX(ICHAN) = 1.0 YMIN(ICHAN) = 0.0 ELSE IF ( NPTS(ICHAN) .EQ. 1 ) THEN DRMS(ICHAN) = 0. DARKRMS(SPTR,ICHAN) = 0. ELSE FLUX(ICHAN) = FLUX(ICHAN) / FLOAT(NPTS(ICHAN)) DARKAVG(SPTR,ICHAN) = DARKAVG(SPTR,ICHAN) $ / FLOAT(NPTS(ICHAN)) DARKOFF(SPTR,ICHAN) = DARKOFF(SPTR,ICHAN) $ / FLOAT(NPTS(ICHAN)) DARKRMS(SPTR,ICHAN) = $ SQRT( DARKRMS(SPTR,ICHAN)/FLOAT(NPTS(ICHAN)) - $ DARKOFF(SPTR,ICHAN)*DARKOFF(SPTR,ICHAN) ) DAVG(ICHAN) = DAVG(ICHAN) / FLOAT(NPTS(ICHAN)) DRMS(ICHAN) = SQRT( DRMS(ICHAN)/FLOAT(NPTS(ICHAN)) $ - DAVG(ICHAN) * DAVG(ICHAN) ) END IF CALL PGSCI ( 5 ) YLO = YPORT(2) - .25*FLOAT(ICHAN) *( YPORT(2)-YPORT(1) ) YHI = YPORT(2) - .25*FLOAT(ICHAN-1)*( YPORT(2)-YPORT(1) ) CALL PGVPORT ( XPORT(1), XPORT(2), YLO, YHI ) WRITE( TITLE,1500) ICHAN, FILTERID(ICHAN) CALL PGMTEXT ( 'L', 3.5, 0.5, 0.5, TITLE ) CALL PGRNGE ( XMIN, XMAX, XLO, XHI ) CALL PGRNGE ( YMIN(ICHAN), YMAX(ICHAN), YLO, YHI ) CALL PGWINDOW ( XLO, XHI, YLO, YHI ) IF ( ICHAN .EQ. 1 ) THEN CALL PGBOX ( 'BCMST', 0.0, 0, 'BCNVST', 0.0, 0 ) ELSE IF ( ICHAN .EQ. 4 ) THEN CALL PGBOX ( 'BCNST', 0.0, 0, 'BCNVST', 0.0, 0 ) ELSE CALL PGBOX ( 'BCST', 0.0, 0, 'BCNVST', 0.0, 0 ) END IF C C Plot the data CALL PGSCI ( 3 ) CALL PGPOINT(N, XPLOT(1), YPLOT(1,ICHAN), SYMBOL(1)) C IF ( N .GT. 1 ) THEN C C And the fit CALL PGSLS ( 2 ) CALL PGLINE (N, XPLOT(1), ZPLOT(1,ICHAN) ) C C And the fit plus the dark offset for this star DO J = 1, N ZPLOT(J,ICHAN) = ZPLOT(J,ICHAN) $ + DARKOFF(SPTR,ICHAN) END DO CALL PGSLS ( 1 ) CALL PGLINE (N, XPLOT(1), ZPLOT(1,ICHAN) ) END IF C---------------------------------------------------------------------- C Edit the data C DAVG(ICHAN) = 0. DRMS(ICHAN) = 0. NPTS(ICHAN) = 0 DO J = 1, N DELTA = ABS( YPLOT(J,ICHAN) - ZPLOT(J,ICHAN) ) C C Scan with missing dark count -- mark the point in red. C IF ( YPLOT(J,ICHAN) .EQ. 0.0 ) THEN C WRITE(9,*) ' EDIT ', J, ' NO DATA' CALL PGSCI(2) CALL PGPOINT(1, XPLOT(J), YPLOT(J,ICHAN), 16 ) C C -- Delete the point if it is in twilight C IF ( ( WHERE(J) .LT. IDUSK ) .OR. $ ( WHERE(J) .GT. IDAWN ) ) THEN C WRITE(9,*) ' TWILIGHT, DELETE POINT ' OKSCAN(WHERE(J)) = $ IOR(OKSCAN(WHERE(J)),MASK(ICHAN) ) C C -- Interpolate good data. Mark the new value in blue. C ELSE C WRITE(9,*) ' INTERPOLATE ' D(WHERE(J),ICHAN) = DARKFIT(WHERE(J),ICHAN) $ + DARKOFF(SPTR, ICHAN) YPLOT(J,ICHAN) = D(WHERE(J),ICHAN) CALL PGSCI(5) CALL PGPOINT(1, XPLOT(J), YPLOT(J,ICHAN), 16) END IF C C Scan with descrepant background -- mark the point in red. C ELSE IF ( DELTA .GT. F*FLUX(ICHAN) ) THEN CALL PGSCI(2) CALL PGPOINT(1, XPLOT(J), YPLOT(J,ICHAN), 16 ) C WRITE(9,*) ' EDIT', J, ' BAD DATA' C C -- Delete points in twilight C IF ( ( WHERE(J) .LT. DUSK ) .OR. $ ( WHERE(J) .GT. DAWN ) ) THEN C WRITE(9,*) ' TWILIGHT, DELETE POINT ' OKSCAN(WHERE(J)) = $ IOR(OKSCAN(WHERE(J)),MASK(ICHAN) ) C C -- Fix the point if we are suppose to. C -- mark the new value with blue. C ELSE IF ( FIXFLAG ) THEN D(WHERE(J),ICHAN) = DARKFIT(WHERE(J),ICHAN) $ + DARKOFF(SPTR, ICHAN) YPLOT(J,ICHAN) = D(WHERE(J),ICHAN) CALL PGSCI(5) CALL PGPOINT(1, XPLOT(J), YPLOT(J,ICHAN), 16) C C -- otherwise delete it. C ELSE C WRITE(9,*) ' DELETE POINT ' OKSCAN(WHERE(J)) = $ IOR(OKSCAN(WHERE(J)),MASK(ICHAN) ) END IF C C Include good scans in the statistics. C ELSE NPTS(ICHAN) = NPTS(ICHAN) + 1 DAVG(ICHAN) = DAVG(ICHAN) + DELTA DRMS(ICHAN) = DRMS(ICHAN) + DELTA * DELTA END IF END DO ! end of editing loop for this channel. C----------------------------------------------------------------------- C C Label right hand side of plot with the dark count measured as a C fraction of the star's brightness. C CALL PGSCI ( 5 ) IF ( FLUX(ICHAN) .GT. 0.0 ) THEN YHI = YHI / FLUX(ICHAN) CALL PGWINDOW ( XLO, XHI, YLO, YHI ) CALL PGBOX ( ' ', 0.0, 0, 'CIMVST', 0.0, 0 ) END IF IF ( NPTS(ICHAN) .EQ. 1 ) THEN DRMS(ICHAN) = 0. ELSE IF ( NPTS(ICHAN) .GT. 1 ) THEN DAVG(ICHAN) = DAVG(ICHAN) / FLOAT( NPTS(ICHAN) ) DRMS(ICHAN) = SQRT ( DRMS(ICHAN)/FLOAT(NPTS(ICHAN)) $ - DAVG(ICHAN) * DAVG(ICHAN) ) END IF END DO ! end of loop over channels CALL PGASK ( .TRUE. ) CALL PGADVANCE 350 CONTINUE CALL PGEND RETURN END SUBROUTINE FIT_2D ( X, Y, YFIT, N, ORDER, RMS, CNT, COEF ) C C Fit a two dimensional polynomial to a data set. C C Input data X(i,j), Y(i), i=1, N, j = 1,2 C Values fitted to the data Y(i), i=1, N C Order of fit ORDER(j), j=1,2 C C RMS = the rms of the fit C COEF the coefficients of the fit C C This subroutine does not output the coefficients of the fit. C C The quality of the fit is evaluated by RMS, the rms of the difference C between the data and the fit. C IMPLICIT UNDEFINED (A-Z) INTEGER *4 N, ORDER(2), CNT REAL *4 X(CNT,2), Y(*), YFIT(*) REAL *4 RMS, MEAN, DELTA INTEGER *4 MAXPAR PARAMETER (MAXPAR=20) REAL *8 A(MAXPAR,MAXPAR), AA(MAXPAR,MAXPAR), WEIGHT, COEF(*) REAL *8 CC(MAXPAR), EE(MAXPAR), BB, VV INTEGER *4 I, J, NUNK SAVE C....................................................................... NUNK = ORDER(1) + ORDER(2) + 1 IF ( ORDER(1) .LT. 0 ) THEN WRITE(6,*) ' ORDER of fit must be non negative ', ORDER(1) ORDER(1) = 0 ELSE IF ( ORDER(2) .LT. 0 ) THEN WRITE(6,*) ' ORDER of fit must be non negative ', ORDER(2) ORDER(2) = 0 ELSE IF ( NUNK+2 .GT. MAXPAR ) THEN WRITE(6,*) ' TOO many parameters. Limit is ',MAXPAR C C Code to change the values of order(i) should go here. C RETURN END IF NUNK = ORDER(1) + ORDER(2) + 1 WEIGHT = 1. CALL LS(1,NUNK,BB,CC,COEF,EE,VV, WEIGHT, A, AA, MAXPAR) MEAN = 0. RMS = 0. DO I = 1, N MEAN = MEAN + Y(I) RMS = RMS + Y(I)*Y(I) BB = Y(I) CC(1) = 1. DO J = 1, ORDER(1) CC(J+1) = X(I,1)**J END DO DO J = 1, ORDER(2) CC(J+1+ORDER(1)) = X(I,2)**J END DO C C LS mode = 2, increment coefficient matrix for least sq fit. C mode = 3, increment solution vector for least sq fit. C The weight is included as real*8 C CALL LS(2,NUNK,BB,CC,COEF,EE,VV, WEIGHT, A, AA, MAXPAR ) CALL LS(3,NUNK,BB,CC,COEF,EE,VV, WEIGHT, A, AA, MAXPAR ) END DO C C LS mode = 4, do the least squares fit C C WRITE(9,*) ' ALL DATA INPUT, READY for FIT...numb of points = ',N DO J = 2, NUNK COEF(J) = 0. END DO IF ( N .EQ. 0 ) THEN COEF(1) = 0. ELSE IF ( N .LT. NUNK ) THEN COEF(1) = MEAN / FLOAT(N) ELSE CALL LS(4,NUNK,BB,CC,COEF,EE,VV, WEIGHT, A, AA, MAXPAR ) END IF MEAN = 0. RMS = 0. DO I = 1, N YFIT(I) = COEF(1) DO J = 1, ORDER(1) YFIT(I) = YFIT(I) + COEF(J+1) * ( X(I,1)**J ) END DO DO J = 1, ORDER(2) YFIT(I) = YFIT(I) + COEF(J+1+ORDER(1)) * ( X(I,2)**J ) END DO DELTA = Y(I) - YFIT(I) MEAN = MEAN + DELTA RMS = RMS + DELTA * DELTA END DO IF ( N .LE. 1 ) THEN RMS = 0. ELSE MEAN = MEAN / FLOAT(N) RMS = SQRT ( RMS / FLOAT(N) - MEAN*MEAN ) END IF RETURN END SUBROUTINE FIT_1D ( X, Y, YFIT, N, ORDER, RMS, COEF ) C C Fit a polynomial to a data set. C C Input data X(i), Y(i), i=1, N C Values fitted to the data Y(i), i=1, N C Order of fit ORDER C Coefficients of the fit COEF C C This subroutine does not output the coefficients of the fit. C C The quality of the fit is evaluated by RMS, the rms of the difference C between the data and the fit. C IMPLICIT UNDEFINED (A-Z) INTEGER *4 N, ORDER, I, J, K, L, NUNK REAL *4 X(*), Y(*), YFIT(*), RMS, MEAN, DELTA INTEGER *4 MAXPAR PARAMETER (MAXPAR=20) REAL *8 A(MAXPAR,MAXPAR), AA(MAXPAR,MAXPAR), WEIGHT, COEF(*) REAL *8 CC(MAXPAR), EE(MAXPAR), BB, VV SAVE C....................................................................... NUNK = ORDER + 1 IF ( NUNK+2 .GT. MAXPAR ) THEN WRITE(6,*) ' TOO many parameters. Limit is ',MAXPAR ORDER = ORDER - 3 NUNK = ORDER + 1 ELSE IF ( ORDER .LT. 0 ) THEN WRITE(6,*) ' ORDER of fit must be greater than zero ', ORDER ORDER = 0 END IF WEIGHT = 1. CALL LS(1,NUNK,BB,CC,COEF,EE,VV, WEIGHT, A, AA, MAXPAR) MEAN = 0. RMS = 0. DO I = 1, N MEAN = MEAN + Y(I) RMS = RMS + Y(I)*Y(I) BB = Y(I) CC(1) = 1.D0 DO J = 1, ORDER CC(J+1) = X(I)**J END DO C WRITE(9,1230) I, X(I), Y(I), (CC(J),J=1,NUNK) C 1230 FORMAT ( ' DATA...', I4, 2F8.3, 2X, 10F8.3 ) C C LS mode = 2, increment coefficient matrix for least sq fit. C mode = 3, increment solution vector for least sq fit. C The weight is included as real*8 C CALL LS(2,NUNK,BB,CC,COEF,EE,VV, WEIGHT, A, AA, MAXPAR ) CALL LS(3,NUNK,BB,CC,COEF,EE,VV, WEIGHT, A, AA, MAXPAR ) END DO C C LS mode = 4, do the least squares fit C C WRITE(9,*) ' ALL DATA INPUT, READY for FIT...numb of points = ',N DO J = 2, NUNK COEF(J) = 0. END DO IF ( N .EQ. 0 ) THEN COEF(1) = 0. ELSE IF ( N .LT. NUNK ) THEN COEF(1) = MEAN / FLOAT(N) ELSE CALL LS(4,NUNK,BB,CC,COEF,EE,VV, WEIGHT, A, AA, MAXPAR ) END IF MEAN = 0. RMS = 0. DO I = 1, N YFIT(I) = COEF(1) DO J = 1, ORDER YFIT(I) = YFIT(I) + COEF(J+1) *( X(I)**J ) END DO DELTA = Y(I) - YFIT(I) MEAN = MEAN + DELTA RMS = RMS + DELTA * DELTA C WRITE(9,1250) I, X(I), Y(I), YFIT(I), DELTA C 1250 FORMAT ( 5X, I5, 3F15.5, 5X, F15.5 ) END DO IF ( N .LE. 1 ) THEN RMS = 0. ELSE MEAN = MEAN / FLOAT(N) RMS = SQRT ( RMS / FLOAT(N) - MEAN*MEAN ) END IF RETURN END SUBROUTINE BOXCAR ( X, Y, YFIT, N, XLO, XHI, YLO, YHI, $ WIDTH, XLINE, YLINE, NLINE, RMS ) C C Do a boxcar smoothing. C C Input data X(i), Y(i), i=1, N C Values fitted to the data Y(i), i=1, N C Range for valid x data XLO, XHI C Range for valid y data YLO, YHI C Width of boxcar in minutes WIDTH C C To facilitate plotting the function, the value of the fit, YLINE(i), C is evaluated on a grid of x values XLINE(i), i=1, NLINE input by C the user. C C The quality of the fit is given by RMS, the rms of the difference C between the data and the fit. C IMPLICIT UNDEFINED (A-Z) INTEGER *4 N, NLINE REAL *4 X(*), Y(*), XLINE(*), YLINE(*), YFIT(*) REAL *4 YLO, YHI, RMS, WIDTH, MEAN, XLO, XHI, DELTA REAL *4 START_X, STOP_X INTEGER *4 I, J, I0, I1, NPTS SAVE C....................................................................... C Convert full width in minutes to half width in hours. C DELTA = WIDTH/120. C....................................................................... C I0 is the first point in range, I1 is the first point beyond range. C I0 = 1 I1 = 1 DO I = 1, N START_X = X(I) - DELTA STOP_X = X(I) + DELTA DO WHILE ( (I0.LT.N) .AND. (X(I0).LT.START_X) ) I0 = I0 + 1 END DO DO WHILE ( (I1.LT.N+1) .AND. (X(I1).LT.STOP_X) ) I1 = I1 + 1 END DO NPTS = I1 - I0 YFIT(I) = 0 IF ( NPTS .GT. 0 ) THEN DO J = I0, I1-1 YFIT(I) = YFIT(I) + Y(J) END DO YFIT(I) = YFIT(I) / NPTS END IF END DO C....................................................................... C Evaluate the boxcar on a grid of points. C I0 = 1 I1 = 1 DO I = 1, NLINE IF ( XLINE(I) .EQ. 0.0 ) THEN YLINE(I) = 0.D0 ELSE START_X = XLINE(I) - DELTA STOP_X = XLINE(I) + DELTA DO WHILE ( (I0.LT.N) .AND. (X(I0).LT.START_X) ) I0 = I0 + 1 END DO DO WHILE ( (I1.LT.N+1) .AND. (X(I1).LT.STOP_X) ) I1 = I1 + 1 END DO NPTS = I1 - I0 YLINE(I) = 0 IF ( NPTS .GT. 0 ) THEN DO J = I0, I1-1 YLINE(I) = YLINE(I) + Y(J) END DO YLINE(I) = YLINE(I) / NPTS ELSE IF ( XLINE(I) .LT. X(1) ) THEN C C Extrapolate the data. Do not interpolate. C YLINE(I) = YFIT(1) ELSE IF ( XLINE(I) .GT. X(N) ) THEN YLINE(I) = YFIT(N) END IF END IF END DO C....................................................................... RETURN END SUBROUTINE LS ( MODE, N, B, C, D, E, VV, WT, A, AA, M ) C C LEAST-SQUARES PROGRAM FOR UP TO 48 VARIABLES WITH C PRECOMPUTED N = NO OF VARIABLES C C changed 18 Nov 86 (DM) maximum number of parameters = 90 C 15 Jan 87 (DM) added weights C 9 Jul 87 (DM) changed defn of VV, so no bomb on PN=N C 5 May 88 (DM) transfered to PC C Eliminated used of unlabeled common. C IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT UNDEFINED (A-Z) INTEGER *4 MODE, N, M, I, J REAL *8 A(M,M), AA(M,M), B, C(M), D(M), E(M), VV REAL *8 WT, PN SAVE C MODE = 1, initialize IF ( N+2 .GT. M ) THEN WRITE(6,*) ' ARRAYS ARE TOO SMALL, LS DIES ' STOP ELSE IF ( MODE .EQ. 1 ) THEN DO 102 I = 1,M DO 101 J = 1,M AA(J,I) = 0. A (J,I) = 0. 101 CONTINUE 102 CONTINUE PN = 0. ELSE IF ( MODE .EQ. 2 ) THEN C MODE = 2, increment coefficients DO 203 I = 1,N DO 202 J = 1,N AA(I,J) = AA(I,J) + C(I)*C(J)*WT 202 CONTINUE 203 CONTINUE ELSE IF ( MODE .EQ. 3 ) THEN C MODE = 3, increment data vector DO 302 I = 1,N A(I,N+1) = A(I,N+1) + B*C(I)*WT A(N+1,I) = A(I,N+1) 302 CONTINUE A(N+1,N+1) = A(N+1,N+1) + B*B * WT PN = PN + 1.D0 ELSE IF ( MODE .EQ. 4 ) THEN C MODE = 4, invert the matrix C A(1,N+2) = 1. A(N+2,1) = -1. DO 410 I = 1,N DO 410 J = 1,N A(I,J) =AA(I,J) 410 CONTINUE CALL INVERT(A, AA, N, M) IF ( PN .LE. DBLE(N) ) THEN VV = 0. ELSE VV = DSQRT(A(1,1)/(PN-DBLE(N))) END IF DO 404 I = 1,N E(I) = VV * SQRT( ABS( A(I+1,I+1) ) ) D(I) = A(I+1,1) 404 CONTINUE ELSE WRITE(6,*) ' LS ERROR: illegal value of MODE = ', MODE END IF RETURN END SUBROUTINE INVERT(A, AA, N, M) C BASICALLY KNOWLES' ROUTINE TO INVERT THE AUGMENTED ARRAY C STORED IN THE LIBRARY NOV 7, 1979 JHS C C changed 18 Nov 86 (DM) max number of parameters = 90 C transfered to PC on 5 May 1988 IMPLICIT UNDEFINED (A-Z) INTEGER *4 MAXPAR, N, M, I, J, K PARAMETER (MAXPAR=20) REAL *8 T(MAXPAR), A(M,M), AA(M,M), Z(MAXPAR,MAXPAR) SAVE IF ( N .GT. M ) THEN WRITE(6,*) ' ARRAYS ARE TOO SMALL. INVERT DIES ' STOP ELSE IF ( M .GT. MAXPAR ) THEN WRITE(6,*) ' MAXPAR IS TOO SMALL, INVERT DIES ' STOP END IF C DO 1 K =1,N DO 2 J = 1,N+2 T(J) = A(1,J) 2 CONTINUE DO 4 I = 1,N+1 DO 3 J = 1,N+1 A(I,J)=A(I+1,J+1)-T(J+1)*A(I+1,1)/T(1) 3 CONTINUE 4 CONTINUE 1 CONTINUE DO 301 J = 2,N+1 DO 300 I = 2,N+1 Z(I,J) = DSQRT(A(I,J)*A(I,J)/(A(I,I)*A(J,J))) 300 CONTINUE 301 CONTINUE RETURN 1004 FORMAT(2X,16F7.3) END