PROGRAM NEWCOMP *- * NEWCOMP * * Reads in a list of spectra, normalises them, and fits a continuum, * and the coadds the, weighting by the inverse variance. * * Input:- List of spectra, including redshifts * * Calls:- RFIG, SVDFIT, PGPLOT * Output:- A file listing the continuum parameters, plots of the fit * * Compile using f77 -C confit.f -lpgplot * * This version reads in errors in the format produced by findnoise.f, * and the continuum you fit is a chebyshev polynomial. Errors * are weighted low in the emission lines and near the ends of the * spectra. * * Header updated 17/8/00, Paul Francis, ANU * Modified to iteratively fit the continuum, modifying an array of * weights as it goes. 22/8/00. *- IMPLICIT NONE INTEGER I, J, K, NBIG INTEGER SCREEN, KEYBD, INFILE, OUTFIL INTEGER COUNT, NPAR, MP, NP INTEGER MAXSIZ, NUM, L, NITER, BIGSIZ PARAMETER (SCREEN=6, BIGSIZ= 20000, : KEYBD=5, INFILE=7, OUTFIL=8, MAXSIZ=5000, MP=5001, NP=10) REAL U(MP,NP), PAR(NP), POL(NP), USIGCLIP, LSIGCLIP REAL WAV(MAXSIZ), SIG(MAXSIZ), ERR(MAXSIZ), INDEX, NORM REAL HIY, Z, V(NP,NP), W(NP), CHISQ, WEIGHT(MAXSIZ) REAL MEAN, MODEL(MAXSIZ), TERR(MAXSIZ) REAL RBWAV(BIGSIZ), RBSIG(BIGSIZ), RBERR(BIGSIZ), VARWT REAL SUMX(BIGSIZ), SUMW(BIGSIZ), NPLOT(BIGSIZ), ZPLOT(BIGSIZ) REAL ZMIN(BIGSIZ), ZMAX(BIGSIZ), WSTART, WEND, WINT REAL CWAV(MAXSIZ), CSIG(MAXSIZ), CERR(MAXSIZ), CN(MAXSIZ) REAL CZ(MAXSIZ), CZMIN(MAXSIZ), CZMAX(MAXSIZ) EXTERNAL POLY CHARACTER*20 ONAME, TNAME, FILNAM * Order of polynomial to fit * WRITE(SCREEN,'(A)') ' Order of polynomial to fit?' * READ(KEYBD,*) NPAR NPAR = 4 * Number of continuum fitting iterations * WRITE(SCREEN,'(A)') ' Number of continuum fitting iterations?' * READ(KEYBD,*) NITER NITER=3 * Upper and lower sigma clipping limits * WRITE(SCREEN,'(A)') ' Upper sigma clipping limit?' * READ(KEYBD,*) USIGCLIP * WRITE(SCREEN,'(A)') ' Lower sigma clipping limit?' * READ(KEYBD,*) LSIGCLIP USIGCLIP = 2.0 LSIGCLIP = 4.0 * Power-law index to fit spectrum to. INDEX = -2.0 * Find the list of spectra: WRITE (SCREEN,'(A)')' List of spectra:' READ (KEYBD,'(A)') FILNAM OPEN (INFILE, FILE=FILNAM, STATUS='OLD') WRITE(SCREEN,'(A)')' ********************************' WRITE(SCREEN,'(A)')' Output file name?' READ(KEYBD,'(A)') ONAME OPEN(OUTFIL, FILE=ONAME, STATUS='NEW') * Set up composite array. WSTART = 1000.0 WINT = 0.5 WEND = 6000.0 NBIG = INT((WEND-WSTART)/WINT)+1 RBWAV(1) = WSTART RBSIG(1) = 0.0 RBERR(1) = 0.0 SUMX(1) = 0.0 SUMW(1) = 0.0 NPLOT(1) = 0.0 ZPLOT(1) = 0.0 ZMIN(1) = 100.0 ZMAX(1) = -1.0 DO 10 I=2, NBIG RBWAV(I) = RBWAV(I-1)+WINT RBSIG(I) = 0.0 RBERR(I) = 0.0 SUMX(I) = 0.0 SUMW(I) = 0.0 NPLOT(I) = 0.0 ZPLOT(I) = 0.0 ZMIN(I) = 100.0 ZMAX(I) = -1.0 10 CONTINUE CALL PGBEGIN(0,'?',1,1) * Loop through all the spectra in the list: DO 110 I=1,1200 * Read in spectrum name: READ(INFILE,20,END=120) TNAME,Z 20 FORMAT(1X,A10,F16.4) * Read in the spectrum IF (Z .GT. 0.0) THEN CALL RFIG(TNAME, NUM, WAV, SIG, ERR, MAXSIZ) * Find the mean flux (and weight the ends of the spectra low) MEAN = 0.0 DO 40 J=1, NUM MEAN = MEAN + SIG(J) IF (WAV(J) .LT. 3400.0 .OR. WAV(J) .GT. 9500.0) THEN ERR(J) = 10.0 ELSE IF (WAV(J) .LT. 3600.0 .OR. WAV(J) .GT. 8500) THEN ERR(J) = 3.0 ELSE IF (WAV(J) .LT. 4000.0 .OR. WAV(J) .GT. 7000) THEN ERR(J) = 2.0 ELSE ERR(J) = 1.0 END IF 40 CONTINUE MEAN = MEAN/REAL(NUM) WRITE(SCREEN,'(1X,I5,1X,A9)') I, TNAME * Scale by the mean and move to rest-frame. Set up initial weights. DO 50 J=1, NUM WAV(J) = WAV(J)/(1.0+Z) SIG(J) = SIG(J)/MEAN WEIGHT(J) = 1.0 IF (WAV(J) .LT. 912.0) THEN WEIGHT(J) = 0.00001 ELSE IF (WAV(J) .GT. 912.0 .AND. WAV(J) .LT. 1180.0) THEN WEIGHT(J) = WEIGHT(J)*0.2 ELSE IF (WAV(J) .GE. 1370.0 .AND. WAV(J) .LT. 1430.0) THEN WEIGHT(J) = WEIGHT(J)*0.001 ELSE IF (WAV(J) .GE. 1480.0 .AND. WAV(J) .LT. 1700.0) THEN WEIGHT(J) = WEIGHT(J)*0.001 ELSE IF (WAV(J) .GE. 1820.0 .AND. WAV(J) .LT. 1950.0) THEN WEIGHT(J) = WEIGHT(J)*0.001 ELSE IF (WAV(J) .GE. 2700.0 .AND. WAV(J) .LT. 2900.0) THEN WEIGHT(J) = WEIGHT(J)*0.001 ELSE IF (WAV(J) .GE. 3700.0 .AND. WAV(J) .LT. 3770.0) THEN WEIGHT(J) = WEIGHT(J)*0.001 ELSE IF (WAV(J) .GE. 4250.0 .AND. WAV(J) .LT. 4450.0) THEN WEIGHT(J) = WEIGHT(J)*0.001 ELSE IF (WAV(J) .GE. 4750.0 .AND. WAV(J) .LT. 5100.0) THEN WEIGHT(J) = WEIGHT(J)*0.001 ELSE IF (WAV(J) .GE. 6300.0 .AND. WAV(J) .LT. 6800.0) THEN WEIGHT(J) = WEIGHT(J)*0.001 END IF 50 CONTINUE * Do the fitting! * Fit a linear combination of the models to the data. Fitting is * done iteratively. DO 80 L=1, NITER DO 55 J=1, NUM TERR(J) = ERR(J)/WEIGHT(J) 55 CONTINUE CALL SVDFIT(WAV,SIG,TERR,NUM,PAR,NPAR,U,V,W,MP,NP, : CHISQ,POLY) * Compute the model, and revise weights. DO 70 J=1,NUM MODEL(J) = 0.0 DO 60 K=1, NPAR CALL POLY(WAV(J),POL,NPAR) MODEL(J) = MODEL(J)+PAR(K)*POL(K) 60 CONTINUE IF (SIG(J) .GT. (MODEL(J)+USIGCLIP*ERR(J)) .OR. : SIG(J) .LT. (MODEL(J)-LSIGCLIP*ERR(J))) THEN WEIGHT(J) = 0.001 ELSE WEIGHT(J) = 1.0 END IF 70 CONTINUE 80 CONTINUE * Divide out the continuum and add a new one in! NORM = 0.0 DO 85 K=1, NPAR CALL POLY(WAV(NUM/2),POL,NPAR) NORM = NORM+PAR(K)*POL(K) 85 CONTINUE NORM = ((WAV(NUM/2)/3000.0)**(-2.0-INDEX))/NORM HIY = -1.0 DO 90 J=1,NUM IF (MODEL(J) .GT. 0.01) THEN SIG(J) = NORM*SIG(J)/MODEL(J) ERR(J) = NORM*ERR(J)/MODEL(J) SIG(J) = SIG(J)*((WAV(J)/3000.0)**(-2.0-INDEX)) ERR(J) = ERR(J)*((WAV(J)/3000.0)**(-2.0-INDEX)) ELSE ERR(J) = 1.0E6 END IF IF (J .GT. (NUM/10)) THEN HIY = MAX(HIY,SIG(J)) END IF 90 CONTINUE * Fill up the composite arrays. COUNT = 1 DO 105 J=1, NBIG IF (RBWAV(J).GT.WAV(1) .AND. RBWAV(J).LT.WAV(NUM-1)) THEN IF (RBWAV(J) .GT. WAV(COUNT+1)) THEN COUNT = COUNT+1 END IF VARWT = 1.0/(ERR(COUNT)**2) SUMX(J) = SUMX(J) + SIG(COUNT)*VARWT SUMW(J) = SUMW(J) + VARWT NPLOT(J) = NPLOT(J) + 1.0 ZPLOT(J) = ZPLOT(J) + Z ZMIN(J) = MIN(ZMIN(J),Z) ZMAX(J) = MAX(ZMAX(J),Z) END IF 105 CONTINUE END IF 110 CONTINUE 120 CLOSE(INFILE) * Assemble the composite HIY = 0.0 DO 200 J=1, NBIG RBSIG(J) = SUMX(J)/SUMW(J) RBERR(J) = 1.0/SUMW(J) ZPLOT(J) = ZPLOT(J)/NPLOT(J) 200 CONTINUE 210 FORMAT(F10.2,2F9.3,I4,3F6.2) * Bin up to something smaller. COUNT = 1 DO 250 J=1,NBIG-5,5 CWAV(COUNT) = 0.0 CSIG(COUNT) = 0.0 CERR(COUNT) = 0.0 CN(COUNT) = NPLOT(J) CZ(COUNT) = ZPLOT(J) CZMIN(COUNT) = ZMIN(J) CZMAX(COUNT) = ZMAX(J) DO 240 K=1, 5 CWAV(COUNT) = CWAV(COUNT) + RBWAV(J+K-1) CSIG(COUNT) = CSIG(COUNT) + RBSIG(J+K-1) CERR(COUNT) = CERR(COUNT) + RBERR(J+K-1) 240 CONTINUE CWAV(COUNT) = CWAV(COUNT)/5.0 CSIG(COUNT) = CSIG(COUNT)/5.0 CERR(COUNT) = CERR(COUNT)/5.0 HIY = MAX(HIY,CSIG(COUNT)) WRITE(OUTFIL,210) CWAV(COUNT), CSIG(COUNT), CERR(COUNT), : INT(CN(COUNT)), CZ(COUNT), CZMIN(COUNT), CZMAX(COUNT) COUNT = COUNT + 1 250 CONTINUE COUNT = COUNT-1 CALL PGSLS(1) CALL PGENV(CWAV(1),CWAV(COUNT),0.0,HIY,0,0) CALL PGLABEL('Rest Wavelength','Flux','Composite') CALL PGLINE(COUNT,CWAV,CSIG) CALL PGSLS(2) CALL PGLINE(COUNT,CWAV,CERR) CLOSE(OUTFIL) CALL PGEND END SUBROUTINE POLY(WAVE,POL,NPAR) *- * Calculates NPARth order polynomial at wavelength WAVE, * and places its various components in the elements of POL * Chebyshev polynomials are used, scaled so that the zero point * is 3000A, with +-1 at 0 and 6000A. *- IMPLICIT NONE INTEGER I, NPAR REAL POL(NPAR), WAVE, TWAV TWAV = (WAVE-3000.0)/3000.0 * Evaluate polynomial iteratively. POL(1) = 1.0 IF (NPAR .GT. 1) THEN POL(2) = TWAV IF (NPAR .GT. 2) THEN DO 10 I=2, NPAR-1 POL(I+1) = 2.0*TWAV*POL(I)-POL(I-1) 10 CONTINUE END IF END IF RETURN END SUBROUTINE RFIG(TNAME, NUM, WAV, SIG, ERRARR, MAXSIZ) *- * RRFIG * Given a character string TNAME, this subroutine reads into the * arrays WAV, SIG and ERR the IUE ascii file data of the spectrum * TNAME.mrg * * Paul Francis, University of Melbourne, Dec 1993 *- INTEGER I, FILNUM, NUM PARAMETER (FILNUM=9) REAL WAV(MAXSIZ), SIG(MAXSIZ), ERRARR(MAXSIZ) REAL WAVE, FLUX CHARACTER TNAME*20, SNAME*21 * Open the file for read access only. SNAME = TNAME OPEN(FILNUM, FILE=SNAME, STATUS='OLD') * Read in the spectrum, discarding zero wavelength rows. NUM = 0 DO 10 I=1,MAXSIZ READ(FILNUM,*,END=20) WAVE, FLUX NUM = NUM + 1 WAV(NUM) = WAVE SIG(NUM) = FLUX ERRARR(NUM) = 1.0 10 CONTINUE * Close the file. 20 CLOSE(FILNUM) END SUBROUTINE SVDFIT(X,Y,SIG,NDATA,A,MA,U,V,W,MP,NP,CHISQ,FUNCS) PARAMETER(NMAX=4001,MMAX=10,TOL=1.0E-5) DIMENSION X(NDATA),Y(NDATA),SIG(NDATA),A(MA),V(NP,NP), : U(MP,NP),W(NP),B(NMAX),AFUNC(NMAX) DO 12 I=1,NDATA CALL FUNCS(X(I),AFUNC,MA) TMP=1.0/SIG(I) DO 11 J=1,MA U(I,J)=AFUNC(J)*TMP 11 CONTINUE B(I)=Y(I)*TMP 12 CONTINUE CALL SVDCMP(U,NDATA,MA,MP,NP,W,V) WMAX=0.0 DO 13 J=1,MA IF(W(J).GT.WMAX)WMAX=W(J) 13 CONTINUE THRESH=TOL*WMAX DO 14 J=1,MA IF(W(J).LT.THRESH)W(J)=0.0 14 CONTINUE CALL SVDKSB(U,W,V,NDATA,MA,MP,NP,B,A) CHISQ=0.0 DO 16 I=1,NDATA CALL FUNCS(X(I),AFUNC,MA) SUM=0.0 DO 15 J=1,MA SUM = SUM+A(J)*AFUNC(J) 15 CONTINUE CHISQ=CHISQ+((Y(I)-SUM)/SIG(I))**2 16 CONTINUE RETURN END SUBROUTINE SVDKSB(U,W,V,M,N,MP,NP,B,X) PARAMETER (NMAX=1000) DIMENSION U(MP,NP),W(NP),V(NP,NP),B(MP),X(NP),TMP(NMAX) DO 12 J=1,N S=0.0 IF (W(J).NE.0.0)THEN DO 11 I=1,M S=S+U(I,J)*B(I) 11 CONTINUE S=S/W(J) ENDIF TMP(J)=S 12 CONTINUE DO 14 J=1,N S=0.0 DO 13 JJ=1,N S=S+V(J,JJ)*TMP(JJ) 13 CONTINUE X(J)=S 14 CONTINUE RETURN END SUBROUTINE SVDCMP(A,M,N,MP,NP,W,V) PARAMETER (NMAX=800) DIMENSION A(MP,NP),W(NP),V(NP,NP),RV1(NMAX) G=0.0 SCALE=0.0 ANORM=0.0 DO 25 I=1,N L=I+1 RV1(I)=SCALE*G G=0.0 S=0.0 SCALE=0.0 IF (I.LE.M) THEN DO 11 K=I,M SCALE=SCALE+ABS(A(K,I)) 11 CONTINUE IF (SCALE.NE.0.0) THEN DO 12 K=I,M A(K,I)=A(K,I)/SCALE S=S+A(K,I)*A(K,I) 12 CONTINUE F=A(I,I) G=-SIGN(SQRT(S),F) H=F*G-S A(I,I)=F-G IF (I.NE.N) THEN DO 15 J=L,N S=0.0 DO 13 K=I,M S=S+A(K,I)*A(K,J) 13 CONTINUE F=S/H DO 14 K=I,M A(K,J)=A(K,J)+F*A(K,I) 14 CONTINUE 15 CONTINUE ENDIF DO 16 K= I,M A(K,I)=SCALE*A(K,I) 16 CONTINUE ENDIF ENDIF W(I)=SCALE *G G=0.0 S=0.0 SCALE=0.0 IF ((I.LE.M).AND.(I.NE.N)) THEN DO 17 K=L,N SCALE=SCALE+ABS(A(I,K)) 17 CONTINUE IF (SCALE.NE.0.0) THEN DO 18 K=L,N A(I,K)=A(I,K)/SCALE S=S+A(I,K)*A(I,K) 18 CONTINUE F=A(I,L) G=-SIGN(SQRT(S),F) H=F*G-S A(I,L)=F-G DO 19 K=L,N RV1(K)=A(I,K)/H 19 CONTINUE IF (I.NE.M) THEN DO 23 J=L,M S=0.0 DO 21 K=L,N S=S+A(J,K)*A(I,K) 21 CONTINUE DO 22 K=L,N A(J,K)=A(J,K)+S*RV1(K) 22 CONTINUE 23 CONTINUE ENDIF DO 24 K=L,N A(I,K)=SCALE*A(I,K) 24 CONTINUE ENDIF ENDIF ANORM=MAX(ANORM,(ABS(W(I))+ABS(RV1(I)))) 25 CONTINUE DO 32 I=N,1,-1 IF (I.LT.N) THEN IF (G.NE.0.0) THEN DO 26 J=L,N V(J,I)=(A(I,J)/A(I,L))/G 26 CONTINUE DO 29 J=L,N S=0.0 DO 27 K=L,N S=S+A(I,K)*V(K,J) 27 CONTINUE DO 28 K=L,N V(K,J)=V(K,J)+S*V(K,I) 28 CONTINUE 29 CONTINUE ENDIF DO 31 J=L,N V(I,J)=0.0 V(J,I)=0.0 31 CONTINUE ENDIF V(I,I)=1.0 G=RV1(I) L=I 32 CONTINUE DO 39 I=N,1,-1 L=I+1 G=W(I) IF (I.LT.N) THEN DO 33 J=L,N A(I,J)=0.0 33 CONTINUE ENDIF IF (G.NE.0.0) THEN G=1.0/G IF (I.NE.N) THEN DO 36 J=L,N S=0.0 DO 34 K=L,M S=S+A(K,I)*A(K,J) 34 CONTINUE F=(S/A(I,I))*G DO 35 K=I,M A(K,J)=A(K,J)+F*A(K,I) 35 CONTINUE 36 CONTINUE ENDIF DO 37 J=I,M A(J,I)=A(J,I)*G 37 CONTINUE ELSE DO 38 J= I,M A(J,I)=0.0 38 CONTINUE ENDIF A(I,I)=A(I,I)+1.0 39 CONTINUE DO 49 K=N,1,-1 DO 48 ITS=1,30 DO 41 L=K,1,-1 NM=L-1 IF ((ABS(RV1(L))+ANORM).EQ.ANORM) GO TO 2 IF ((ABS(W(NM))+ANORM).EQ.ANORM) GO TO 1 41 CONTINUE 1 C=0.0 S=1.0 DO 43 I=L,K F=S*RV1(I) IF ((ABS(F)+ANORM).NE.ANORM) THEN G=W(I) H=SQRT(F*F+G*G) W(I)=H H=1.0/H C= (G*H) S=-(F*H) DO 42 J=1,M Y=A(J,NM) Z=A(J,I) A(J,NM)=(Y*C)+(Z*S) A(J,I)=-(Y*S)+(Z*C) 42 CONTINUE ENDIF 43 CONTINUE 2 Z=W(K) IF (L.EQ.K) THEN IF (Z.LT.0.0) THEN W(K)=-Z DO 44 J=1,N V(J,K)=-V(J,K) 44 CONTINUE ENDIF GO TO 3 ENDIF IF (ITS.EQ.30) PAUSE 'No convergence in 30 iterations' X=W(L) NM=K-1 Y=W(NM) G=RV1(NM) H=RV1(K) F=((Y-Z)*(Y+Z)+(G-H)*(G+H))/(2.0*H*Y) G=SQRT(F*F+1.0) F=((X-Z)*(X+Z)+H*((Y/(F+SIGN(G,F)))-H))/X C=1.0 S=1.0 DO 47 J=L,NM I=J+1 G=RV1(I) Y=W(I) H=S*G G=C*G Z=SQRT(F*F+H*H) RV1(J)=Z C=F/Z S=H/Z F= (X*C)+(G*S) G=-(X*S)+(G*C) H=Y*S Y=Y*C DO 45 NM=1,N X=V(NM,J) Z=V(NM,I) V(NM,J)= (X*C)+(Z*S) V(NM,I)=-(X*S)+(Z*C) 45 CONTINUE Z=SQRT(F*F+H*H) W(J)=Z IF (Z.NE.0.0) THEN Z=1.0/Z C=F*Z S=H*Z ENDIF F= (C*G)+(S*Y) X=-(S*G)+(C*Y) DO 46 NM=1,M Y=A(NM,J) Z=A(NM,I) A(NM,J)= (Y*C)+(Z*S) A(NM,I)=-(Y*S)+(Z*C) 46 CONTINUE 47 CONTINUE RV1(L)=0.0 RV1(K)=F W(K)=X 48 CONTINUE 3 CONTINUE 49 CONTINUE RETURN END