program Tempfit c c program reset to read only one template from fort.28 (vel, f_lambda/f_max) c + read all fitted profile from fort.27 (vel, f_lambda), and only points c relevant for the fit from fort.26 (vel, f_lambda, err_lambda) c c This program makes a best chi**2 fit to a blend using a template. The template c is the blue wing of Ly alpha, read from fort.28, which has c the profiles from the best fit Gaussians model (+ the data normalized to 1 in c flux, wavelength is in velocity units). The data points which are to be taken c into account when calculating the chi**2 are read from fort.26. c The fit parameters are read from fort.33. c Output: fort.2- best fit model + all the data points c fort.4- data points used in the fit, delta chi. c To fit a new object: 1. load ly alpha template to fort.28 c To fit a new line: 1.load lines to fort.26 + fort.27 2. reset fort.33 c check below what line ratio is asummed (iequal in MRQMIN) c Note- best fit solution can be very sensitive to intial guess c and different compilers sometimes give different results(!) c parameter (ne=2000,ne1=4000,me=30) integer lista(me) real fline(ne,3),a(me),covar(me,me),alpha(me,me),xwav(ne) & ,vel1(ne),flz1(ne),sig1(ne),tmpl(ne,2),yt(ne) & , yc1(ne),yc2(ne),yc3(ne),yc4(ne),yc5(ne),ycomp(ne,10) common /fga/ifix,deltav,deltav1,iequal common /tmp/tmpl,ltemp external ftemp iequal=0 if (iequal.eq.0) print*,'assume optically thin case (2:1)' if (iequal.eq.1) print*,'assume optically thick case (1:1)' c c iequal= 1- use equal strength in doublets, iequal= 0- use 2:1 c ifix= 1 - assumes first component is a doublet. Use for O VI, C IV, Mg II. c ifix= 2 - assumes 2nd +3rd components are a doublet. Use for Ly alpha c ifix= 3 - assumes 1+2 and 3+4 are both doublets. Can be used for c C IV + Si II 1526 1533 c ifix= 5 - used for the Si IV doublet + O IV] quintuplet c read(33,*) wline, ifix, ma, mfit read(33,*) (lista(i),i=1,ma) if(ifix.ne.5) then read(33,*) (a(i),i=1,ma) else read(33,*) a(1),a(7),a(2) a(3)=a(1)/2. a(5)=a(7)*2 a(9)=a(7)*6 a(11)=a(7)*4 a(13)=a(7)*2 a(4)=a(2)+1932.3 a(8)=a(2)+1291.4 a(6)=a(8)-546.1 a(10)=a(8)+295.2 a(12)=a(8)+1077.3 a(14)=a(8)+1629. endif if(ifix.ne.3) read(33,*) deltav if(ifix.eq.3) read(33,*) deltav,deltav1 write(6,*) wline, ifix, ma, mfit write(6,*) (lista(i),i=1,ma) c c a(1)-line peak flux density, a(2)- line velocity shift c write(6,'(1p,E11.3,0p,f11.2)') (a(i),i=1,ma) if(ifix.ne.3) write(6,*) deltav if(ifix.eq.3) write(6,*) deltav,deltav1 c-- read data cleaned of absorption, restrict fitted range to v1-v2 read(33,*) v1, v2 read(33,*) z write(6,*) 'z=',z c-- make fit to the Si IV + O IV] 1400A multiplet (2+5 components) iv=0 j=0 c c read all data points of the fitted line profile c do 1 i=1,10000 1 read(27,*,end=109) fline(i,1),fline(i,2) c c read data points of the fitted line profile which are used in the fit c 109 lline=i-1 do 11 i=1,10000 read(26,*,end=110) x1,x2,x3 if(x1.lt.v1.or.x1.gt.v2) goto 11 iv=iv+1 vel1(iv)=x1 flz1(iv)=x2 sig1(iv)=x3 11 continue 110 ndd=iv print*,ndd,v1,v2 c c read template profile c do 2 i=1,10000 read(28,*,end=111) tmpl(i,1),tmpl(i,2) 2 continue 111 ltemp=i-1 k=0 c c making tmpl symmetric by using the blue wing only c do 3 i=1,ltemp c-- do not change negative velocity part if(tmpl(i,1).le.0.) goto 3 k=k+1 c-- im= first bin of the red wing k= bin count in red wing only if(k.eq.1) im=i c-- if red wing is longer than blue wing, end wing. if(im-k.le.0) goto 33 tmpl(i,1)=-tmpl(im-k,1) tmpl(i,2)=tmpl(im-k,2) 3 continue 33 ltemp=i-1 ftotal=0. c c calculating integrated area of template (converted later to erg/cm^2/sec) c do 9 i=2,ltemp-1 dv=(tmpl(i+1,1)-tmpl(i-1,1))/2. df=dv*tmpl(i,2) ftotal=ftotal+df 9 continue alamda=-5. chi=0. nad=5 44 do 4 i=1,nad call mrqmin(vel1,flz1,sig1,ndd,a,ma,lista,mfit,covar,alpha, & mfit,chisq,ftemp,alamda) if(i.eq.nad-2.or.i.eq.nad) print*,i,chisq 4 continue print*,'type in # of additional iterations' read*,nad if(nad.gt.0) goto 44 ytmax=0. do 5 i=1,lline yt(i)=0. c c xwav - units of wavelength c xwav(i)=(fline(i,1)/3.e5+1.)*wline do 6 j=1,ma-1,2 ncomp=(j+1)/2 ycomp(i,ncomp)=a(j)*tmpflx(fline(i,1)-a(j+1)) yt(i)=yt(i)+ycomp(i,ncomp) if(yt(i).gt.ytmax) ytmax=yt(i) 6 continue if(ifix.eq.1) then yc1(i)=ycomp(i,1)+ycomp(i,2) yc2(i)=ycomp(i,3) yc3(i)=ycomp(i,4) yc4(i)=ycomp(i,5) yc5(i)=ycomp(i,6) endif if(ifix.eq.2) then yc1(i)=ycomp(i,1) yc2(i)=ycomp(i,2)+ycomp(i,3) endif if(ifix.eq.3) then yc1(i)=ycomp(i,1)+ycomp(i,2) yc2(i)=ycomp(i,3) yc3(i)=ycomp(i,4) endif if(ifix.eq.5) then yc1(i)=ycomp(i,1)+ycomp(i,2) yc2(i)=ycomp(i,3)+ycomp(i,4)+ycomp(i,5)+ycomp(i,6)+ & ycomp(i,7) endif 5 continue c c print observed and fit profiles c do 12 i=1,lline c c find FWHM of model fit c yt(i)=yt(i)/ytmax if(i.eq.1) goto 12 ytemp1=tmpflx(fline(i-1,1)) ytemp2=tmpflx(fline(i,1)) sign=(yt(i)-0.5)*(yt(i-1)-0.5) if(sign.le.0.) then aa=(fline(i,1)-fline(i-1,1))/(yt(i)-yt(i-1)) bb=fline(i,1)-aa*yt(i) if(vhalf1.eq.0.) vhalf1=aa*0.5+bb if(vhalf1.ne.0.) vhalf2=aa*0.5+bb endif sign=(ytemp2-0.5)*(ytemp1-0.5) if(sign.le.0.) then aa=(fline(i,1)-fline(i-1,1))/(ytemp2-ytemp1) bb=fline(i,1)-aa*ytemp2 if(vtemp1.eq.0.) vtemp1=aa*0.5+bb if(vtemp1.ne.0.) vtemp2=aa*0.5+bb endif if(ifix.eq.1) write(2,'(f9.2,7g12.4)')xwav(i),fline(i,2), & yt(i)*ytmax,yc1(i),yc2(i),yc3(i),yc4(i),yc5(i) if(ifix.eq.2) write(2,'(f9.2,5g12.4)')xwav(i),fline(i,2), & yt(i)*ytmax,yc1(i),yc2(i),ycomp(i,4) if(ifix.eq.3) write(2,'(f9.2,5g12.4)')xwav(i),fline(i,2), & yt(i)*ytmax,yc1(i),yc2(i),yc3(i) if(ifix.eq.5) write(2,'(f9.2,5g12.4)')xwav(i),fline(i,2), & yt(i)*ytmax,yc1(i),yc2(i) cc write(2,'(f11.2,2f8.4)') fline(i,1),fline(i,2)/ytmax,yt(i) 12 continue c c print observed profile and deviation from fit in units of sigma c do 7 i=1,ndd yt(i)=0. do 8 j=1,ma-1,2 yt(i)=yt(i)+a(j)*tmpflx(vel1(i)-a(j+1)) 8 continue cc write(4,*) vel1(i),flz1(i),(flz1(i)-yt(i))/sig1(i) 7 continue print*,'FWHM=',vhalf2-vhalf1 print*,'FWHM Ly alpha=',vtemp2-vtemp1 c c total line flux ftotline = wline*(1+z.)*a(i)/c, where i=1,3,5.. c print each components v, total flux c do 10 i=1,ma,2 ftotline=ftotal*a(i)*wline*(1.+z)/3.e5 print*, ftotline, a(i+1) 10 continue end SUBROUTINE ftemp(X,A,Y,DYDA,NA) c c sums up na/2 symmetric Ly alpha profiles c DIMENSION A(NA),DYDA(NA) Y=0. DO 1 I=1,NA-1,2 ex=tmpflx(x-a(i+1)) ex1=tmpflx(x-(a(i+1)-3.)) ex2=tmpflx(x-(a(i+1)+3.)) y=y+a(i)*ex DYDA(I)=ex DYDA(I+1)=a(i)*(ex2-ex1)/6. 1 CONTINUE cc write(10,*) x,y RETURN END function tmpflx(x) parameter (ne=2000,ne1=4000,me=30) real tmpl(ne,2) common /tmp/tmpl,ltemp c c-- if x (=velocity) is outside range then flux=0. c tmpflx=0. do 1 j=1,ltemp-1 if(tmpl(j,1).lt.x.and.tmpl(j+1,1).ge.x)then d=(x-tmpl(j,1))/(tmpl(j+1,1)-tmpl(j,1)) tmpflx=(1.-d)*tmpl(j,2)+d*tmpl(j+1,2) endif 1 continue return end SUBROUTINE MRQMIN(X,Y,SIG,NDATA,A,MA,LISTA,MFIT, * COVAR,ALPHA,NCA,CHISQ,FUNCS,ALAMDA) PARAMETER (MMAX=30) DIMENSION X(NDATA),Y(NDATA),SIG(NDATA),A(MA),LISTA(MA), * COVAR(NCA,NCA),ALPHA(NCA,NCA),ATRY(MMAX),BETA(MMAX),DA(MMAX) common /fga/ifix,deltav,deltav1,iequal IF(ALAMDA.LT.0.)THEN KK=MFIT+1 DO 12 J=1,MA IHIT=0 DO 11 K=1,MFIT IF(LISTA(K).EQ.J)IHIT=IHIT+1 11 CONTINUE IF (IHIT.EQ.0) THEN LISTA(KK)=J KK=KK+1 ELSE IF (IHIT.GT.1) THEN PAUSE 'Improper permutation in LISTA' ENDIF 12 CONTINUE IF (KK.NE.(MA+1)) PAUSE 'Improper permutation in LISTA' ALAMDA=0.001 CALL MRQCOF(X,Y,SIG,NDATA,A,MA,LISTA,MFIT,ALPHA,BETA,NCA,CHISQ,F *UNCS) OCHISQ=CHISQ DO 13 J=1,MA ATRY(J)=A(J) 13 CONTINUE ENDIF DO 15 J=1,MFIT DO 14 K=1,MFIT COVAR(J,K)=ALPHA(J,K) 14 CONTINUE COVAR(J,J)=ALPHA(J,J)*(1.+ALAMDA) DA(J)=BETA(J) 15 CONTINUE CALL GAUSSJ(COVAR,MFIT,NCA,DA,1,1) c-- option for alamda=0 was deleted (covsrt is not used) DO 16 J=1,MFIT ATRY(LISTA(J))=A(LISTA(J))+DA(J) c-- following condition prevents negative fluxes cc if(j.gt.6) goto 16 j2=j/2 if(j2*2.ne.j) goto 16 if(ATRY(LISTA(J)).lt.0.) ATRY(LISTA(J))=ATRY(LISTA(J))/100. 16 CONTINUE if(iequal.eq.1) then if(ifix.eq.1.or.ifix.eq.3) atry(3)=atry(1) if(ifix.eq.3) atry(7)=atry(5) if(ifix.eq.2) atry(5)=atry(3) else if(ifix.eq.1.or.ifix.eq.3) atry(3)=atry(1)/2. if(ifix.eq.3) atry(7)=atry(5)/2. if(ifix.eq.2) atry(5)=atry(3)/2. endif if(ifix.eq.1.or.ifix.eq.3) atry(4)=atry(2)+deltav if(ifix.eq.3) atry(8)=atry(6)+deltav1 if(ifix.eq.2) atry(6)=atry(4)+deltav if(ifix.eq.5) then atry(3)=atry(1)/2. atry(5)=atry(7)*2 atry(9)=atry(7)*6 atry(11)=atry(7)*4 atry(13)=atry(7)*2 atry(4)=atry(2)+1932.3 atry(8)=atry(2)+1291.4 atry(6)=atry(8)-546.1 atry(10)=atry(8)+295.2 atry(12)=atry(8)+1077.3 atry(14)=atry(8)+1629. endif cc write(6,*) (atry(i),i=1,14) CALL MRQCOF(X,Y,SIG,NDATA,ATRY,MA,LISTA,MFIT,COVAR,DA,NCA,CHISQ,FU *NCS) IF(CHISQ.LT.OCHISQ)THEN ALAMDA=0.1*ALAMDA OCHISQ=CHISQ DO 18 J=1,MFIT DO 17 K=1,MFIT ALPHA(J,K)=COVAR(J,K) 17 CONTINUE BETA(J)=DA(J) A(LISTA(J))=ATRY(LISTA(J)) cc if(j.gt.6) goto 18 j2=j/2 if(j2*2.ne.j) goto 18 if(A(LISTA(J)).lt.0.) A(LISTA(J))=A(LISTA(J))/100. 18 CONTINUE if(iequal.eq.1) then if(ifix.eq.1.or.ifix.eq.3) a(3)=a(1) if(ifix.eq.3) a(7)=a(5) if(ifix.eq.2) a(5)=a(3) else if(ifix.eq.1.or.ifix.eq.3) a(3)=a(1)/2. if(ifix.eq.3) a(7)=a(5)/2. if(ifix.eq.2) a(5)=a(3)/2. endif if(ifix.eq.1.or.ifix.eq.3) a(4)=a(2)+deltav if(ifix.eq.3) a(8)=a(6)+deltav1 if(ifix.eq.2) a(6)=a(4)+deltav if(ifix.eq.5) then a(3)=a(1)/2. a(5)=a(7)*2 a(9)=a(7)*6 a(11)=a(7)*4 a(13)=a(7)*2 a(4)=a(2)+1932.3 a(8)=a(2)+1291.4 a(6)=a(8)-546.1 a(10)=a(8)+295.2 a(12)=a(8)+1077.3 a(14)=a(8)+1629. endif ELSE ALAMDA=10.*ALAMDA CHISQ=OCHISQ ENDIF RETURN END SUBROUTINE MRQCOF(X,Y,SIG,NDATA,A,MA,LISTA,MFIT,ALPHA,BETA,NALP,CH *ISQ,FUNCS) PARAMETER (MMAX=30) external funcs DIMENSION X(NDATA),Y(NDATA),SIG(NDATA),ALPHA(NALP,NALP),BETA(MA), * DYDA(MMAX),LISTA(MFIT),A(MA) DO 12 J=1,MFIT DO 11 K=1,J ALPHA(J,K)=0. 11 CONTINUE BETA(J)=0. 12 CONTINUE CHISQ=0. DO 15 I=1,NDATA CALL FUNCS(X(I),A,YMOD,DYDA,MA) SIG2I=1./(SIG(I)*SIG(I)) DY=Y(I)-YMOD DO 14 J=1,MFIT WT=DYDA(LISTA(J))*SIG2I DO 13 K=1,J ALPHA(J,K)=ALPHA(J,K)+WT*DYDA(LISTA(K)) 13 CONTINUE BETA(J)=BETA(J)+DY*WT 14 CONTINUE CHISQ=CHISQ+DY*DY*SIG2I 15 CONTINUE DO 17 J=2,MFIT DO 16 K=1,J-1 ALPHA(K,J)=ALPHA(J,K) 16 CONTINUE 17 CONTINUE RETURN END SUBROUTINE GAUSSJ(A,N,NP,B,M,MP) PARAMETER (NMAX=50) DIMENSION A(NP,NP),B(NP,MP),IPIV(NMAX),INDXR(NMAX),INDXC(NMAX) DO 11 J=1,N IPIV(J)=0 11 CONTINUE DO 22 I=1,N BIG=0. DO 13 J=1,N IF(IPIV(J).NE.1)THEN DO 12 K=1,N IF (IPIV(K).EQ.0) THEN IF (ABS(A(J,K)).GE.BIG)THEN BIG=ABS(A(J,K)) IROW=J ICOL=K ENDIF ELSE IF (IPIV(K).GT.1) THEN PAUSE 'Singular matrix' ENDIF 12 CONTINUE ENDIF 13 CONTINUE IPIV(ICOL)=IPIV(ICOL)+1 IF (IROW.NE.ICOL) THEN DO 14 L=1,N DUM=A(IROW,L) A(IROW,L)=A(ICOL,L) A(ICOL,L)=DUM 14 CONTINUE DO 15 L=1,M DUM=B(IROW,L) B(IROW,L)=B(ICOL,L) B(ICOL,L)=DUM 15 CONTINUE ENDIF INDXR(I)=IROW INDXC(I)=ICOL IF (A(ICOL,ICOL).EQ.0.) PAUSE 'Singular matrix.' PIVINV=1./A(ICOL,ICOL) A(ICOL,ICOL)=1. DO 16 L=1,N A(ICOL,L)=A(ICOL,L)*PIVINV 16 CONTINUE DO 17 L=1,M B(ICOL,L)=B(ICOL,L)*PIVINV 17 CONTINUE DO 21 LL=1,N IF(LL.NE.ICOL)THEN DUM=A(LL,ICOL) A(LL,ICOL)=0. DO 18 L=1,N A(LL,L)=A(LL,L)-A(ICOL,L)*DUM 18 CONTINUE DO 19 L=1,M B(LL,L)=B(LL,L)-B(ICOL,L)*DUM 19 CONTINUE ENDIF 21 CONTINUE 22 CONTINUE DO 24 L=N,1,-1 IF(INDXR(L).NE.INDXC(L))THEN DO 23 K=1,N DUM=A(K,INDXR(L)) A(K,INDXR(L))=A(K,INDXC(L)) A(K,INDXC(L))=DUM 23 CONTINUE ENDIF 24 CONTINUE RETURN END