character*300 incat,icomm*500,ofil*400 character*20 icol(20),chs*5 integer ico(20) double precision p(200000,100),par(1000000),acc, +snlim,mlim,b,c,mcut,flux real xv1(200000),xv2(200000),xmi,xma,xb(200000), +yb(200000),eyb(200000),yb2(200000),xbfit(200000), +ybfit(200000), +xi(200000),yi(200000),xx(200000),yy(200000) integer ip(200000) external bin_mat,FCHI,CSPIN common /CSTAK/DSTAK(1000000) call ISTKIN(1000000,4) ! reading catalog filename read *, incat ! defining column names ncol=12 icol(1)='NUMBER' icol(2)='X_IMAGE' icol(3)='Y_IMAGE' icol(4)='KRON_RADIUS' icol(5)='MAG_AUTO' icol(6)='MAGERR_AUTO' icol(7)='FWHM_IMAGE' icol(8)='CLASS_STAR' icol(9)='FLUX_APER' icol(10)='FLUXERR_APER' icol(11)='IMAFLAGS_ISO' icol(12)='FLAGS' ! copy header into separate file icomm='grep \# '//incat(1:len_trim(incat))// +' > temp_header' call system(icomm) ! associate col number to col name do i=1,ncol icomm='grep '//icol(i)(1:len_trim(icol(i))) +//' temp_header '// +' | awk ''{print $2,$3}'' > temp1' call system(icomm) open(21,file='temp1',status='old') read(21,*) ico(i) close(21) enddo iid=ico(1) ixce=ico(2) iyce=ico(3) irkron=ico(4) imkron=ico(5) iemkron=ico(6) ifwhm=ico(7) isi=ico(8) imag=ico(9)+4 iemag=ico(10)+4 iflagi=ico(11) iflag=ico(12) ! number of columns to read ic=max(iid,ixce,iyce,irkron,imkron,iemkron,ifwhm,isi,imag, +iemag,iflagi,iflag) ! counting number of header lines icomm='wc -l temp_header > temp1' call system(icomm) open(21,file='temp1',status='old') read(21,*) nheader close(21) ! reading the catalog open(21,file=incat,status='old') do i=1,nheader read(21,*) enddo i=0 100 i=i+1 read(21,*,END=110) (p(i,j),j=1,ic) goto 100 110 n=i-1 close(21) c defining vectors with S/N and 2arcsec diam. mag. nx=0 do i=1,n if ((p(i,imag).gt.0.D+00).and.(NINT(p(i,iflagi)).eq.0) +.and.(NINT(p(i,iflag)).eq.0).and. +(( -2.5D+00*DLog10(p(i,imag)) ).lt.30.D+00)) then c if (p(i,isi).gt.0.8D+00) then nx=nx+1 xv1(nx)=sngl( -2.5D+00*DLog10(p(i,imag)) ) xv2(nx)=sngl( p(i,imag)/p(i,iemag) ) c endif endif enddo c binning the S/N vector according to mag, each bin including the same c number of objects np, with min allowed number=npmin np=30 npmin=10 call bin_mat(nx,xv1,xv2,np,npmin,nb,xb,yb,eyb,yb2,ip) c fitting the binned trend with the analytic expression c p1*flux/(p2*flux+p3) mlim=0.D+00 304 acc=1.D-05 nmax=10000 if (mlim.eq.0.D+00) then mcut=22.D+00 else mcut=mlim-4.D+00 endif do i=1,nb c only faint objects (>=21) are used in the fit if (xb(i).gt.sngl(mcut)) then par(200000+i)=1.D+00 else par(200000+i)=0.D+00 endif par(10+i)=dble(xb(i)) par(100000+i)=dble(yb(i)) enddo par(10)=DFLOAT(nb) if (mlim.eq.0.D+00) then c computing initial conditions for p1, p2, and p3 par(1)=1.D+00 par(2)=10.D+00**(-0.4D+00*(par(10+1)-25.D+00))/par(100000+1)**2 par(3)=par(1)*10.D+00**(-0.4D+00*(par(10+nb)-25.D+00))/ +par(100000+nb) par(3)=par(3)**2 endif npar=10 call DSMNF(npar,par,FCHI,nmax,acc) snlim=5.D+00 b=-(snlim/par(1))**2*par(2) c=-(snlim/par(1))**2*par(3) if (mlim.eq.0.D+00) then c repeat the minimization using a fixed mag range mlim=(-b+dsqrt(b**2-4.D+00*c))/2.D+00 mlim=-2.5D+00*DLog10(mlim)+25.D+00 goto 304 endif mlim=(-b+dsqrt(b**2-4.D+00*c))/2.D+00 mlim=-2.5D+00*DLog10(mlim)+25.D+00 ccc if binned trend goes deep enough, we use it to compute mlim ccc if (yb(nb).lt.5.D+00) then ccc call SRTPAR(yb,1,ip,1,nb) ccc ii=0 ccc do i=1,nb ccc if (xb(ip(i)).ge.(mlim-2.)) then ccc ii=ii+1 ccc xx(ii)=yb(ip(i)) ccc yy(ii)=xb(ip(i)) ccc endif ccc enddo ccc xi(1)=5.E+00 ccc call CSPIN(xx,yy,ii,xi,yi,1) ccc mlim=yi(1) ccc endif c set limits for the plot xmi=mlim-3.5 xma=mlim+1. ymi=0.1 yma=99. c compute fitted function k=0 505 k=k+1 xbfit(k)=xmi+0.02*FLOAT(k-1) flux=10.D+00**(-0.4D+00*(dble(xbfit(k))-25.D+00)) ybfit(k)=sngl(par(1)*flux/DSqRt(flux*par(2)+par(3))) if (xbfit(k).lt.xma) goto 505 nfit=k c finally we do the plot ofil='MLIM.ps' call PGOPEN('PLOTS/'//ofil(1:len_trim(ofil))// +'/VCPS') c call PGOPEN('1/XS') call PGPAP(0,1.) call PGSCR(0,1.,1.,1.) call PGSCR(1,0.,0.,0.) call PGSLW(2) call PGSCH(1.2) call PGENV(xmi,xma,ymi,yma,0,0) call PGLAB('mag (2'''' diameter)','S/N','') call PGSCH(0.45) call PGLAB('','',incat(1:(len_trim(incat)-4))) call PGSCH(1.) call PGPT(nx,xv1,xv2,1) call PGSCI(2) call PGLINE(nb,xb,yb) call PGTEXT(xma-0.6*(xma-xmi),yma-0.1*(yma-ymi),'binned trend') call PGSLW(3) call PGSCI(4) call PGLINE(nfit,xbfit,ybfit) call PGSLW(2) call PGTEXT(xma-0.6*(xma-xmi),yma-0.17*(yma-ymi),'best-fit') c call PGLINE(nb,xb,yb2) call PGSCI(1) write(chs,FMT='(F5.2)') mlim call PGTEXT(xma-0.6*(xma-xmi),yma-0.24*(yma-ymi), +'mag\dlim\u='//chs) call PGSLS(2) call PGMOVE(xmi,sngl(snlim)) call PGDRAW(xma,sngl(snlim)) call PGSLS(1) call PGMOVE(sngl(mlim),ymi) call PGDRAW(sngl(mlim),yma) call PGEND c store mag_lim and fitting parameters in MLIM.dat open(21,FILE='OUTPUT_RES/MLIM.dat',status='unknown') write(21,FMT='(A30)') '#SN_LIM MLIM p_i ' snlim=5.D+00 b=-(snlim/par(1))**2*par(2) c=-(snlim/par(1))**2*par(3) mlim=(-b+dsqrt(b**2-4.D+00*c))/2.D+00 mlim=-2.5D+00*DLog10(mlim)+25.D+00 write(21,FMT='(F6.1,F8.3,E16.6)') snlim,mlim,par(1) snlim=10.D+00 b=-(snlim/par(1))**2*par(2) c=-(snlim/par(1))**2*par(3) mlim=(-b+dsqrt(b**2-4.D+00*c))/2.D+00 mlim=-2.5D+00*DLog10(mlim)+25.D+00 write(21,FMT='(F6.1,F8.3,E16.6)') snlim,mlim,par(2) snlim=15.D+00 b=-(snlim/par(1))**2*par(2) c=-(snlim/par(1))**2*par(3) mlim=(-b+dsqrt(b**2-4.D+00*c))/2.D+00 mlim=-2.5D+00*DLog10(mlim)+25.D+00 write(21,FMT='(F6.1,F8.3,E16.6)') snlim,mlim,par(3) close(21) c ofil='PLOTS/M'//incat(1:(len_trim(incat)-4))//'.ps' c icomm='mv -f MLIM.ps '//ofil(1:len_trim(ofil)) c call system(icomm) stop end c subroutine FCHI(npar,par,nf,f) double precision par(*),f,flux nb=NINT(par(10)) f=0.D+00 do i=1,nb flux=10.D+00**(-0.4D+00*(par(10+i)-25.D+00)) par(300000+i)=par(1)*flux/DSqRt(flux*par(2)+par(3)) f=f+par(200000+i)*(par(100000+i)-par(300000+i))**2 enddo f=f/DFLOAT(nb) return end c c c subroutine bin_mat(n,xv1,xv2,np,npmin,nb,xb,yb,eyb,yb2,ip) real xv1(*),xv2(*),xb(*),yb(*),eyb(*),yb2(*) integer ip(*) double precision v(200000),cc,t,s external bwsm,median,SRTPAR call SRTPAR(xv1,1,ip,1,n) i1=-np+1 i=0 643 i=i+1 i1=i1+np i2=i1+np-1 if ((i1.ge.n).or.((n-i1).lt.npmin)) then goto 646 endif i2=min(i2,n) x1=xv1(ip(i1)) x2=xv1(ip(i2)) xb(i)=(x1+x2)/2. c write(6,*) i1,i2,n,xb(i) do ii=i1,i2 v(ii-i1+1)=dble( xv2(ip(ii)) ) c if (xb(i).gt.23.6) then c write(6,*) xv1(ip(ii)),xv2(ip(ii)),xb(i) c endif enddo cc=1.5D+00 call bwsm(i2-i1+1,v,cc,t,s) yb(i)=sngl(t) eyb(i)=sngl(s) c write(6,*) xb(i),yb(i) call median(i2-i1+1,v,t) yb2(i)=sngl(t) if (i1.lt.n) goto 643 646 nb=i-1 return end c c subroutine fai_stat(kk,xx,cc,theta,etheta,sigma,esigma,med,emed $ ,med2,emed2,sigma2,esigma2,skew,eskew1,eskew2,nsim) double precision xx(*),cc,theta,etheta,sigma,esigma,med,emed,c,t,s double precision vv(10000000),vt(10000),vs(10000),vm(10000),med2 $ ,sigma2,emed2,esigma2,vt2(10000),vs2(10000),skew,eskew1 $ ,eskew2 ,vsk(10000) external bwsm,medstd2 do i=1,nsim do k=1,kk if (i.eq.1) then kv=k else kv=IGNUIN(1,kk) endif vv(k)=xx(kv) enddo call stat(kk,vv,skew) c=cc t=-1.D+00 s=-1.D+00 call bwsm(kk,vv,c,t,s) call median(kk,vv,med) call medstd2(kk,vv,med2,sigma2) vt(i)=t vs(i)=s vm(i)=med vt2(i)=med2 vs2(i)=sigma2 vsk(i)=skew c write(6,*) skew c read(5,*) enddo c write(6,*) vt(1),c c stop c=2.5D+00 t=-1.D+00 s=-1.D+00 call bwsm(nsim,vt,c,t,s) theta=t etheta=s t=-1.D+00 s=-1.D+00 call bwsm(nsim,vs,c,t,s) sigma=vs(1) esigma=s/t*vs(1) t=-1.D+00 s=-1.D+00 call bwsm(nsim,vm,c,t,s) med=t emed=s t=-1.D+00 s=-1.D+00 call bwsm(nsim,vt2,c,t,s) med2=t emed2=s t=-1.D+00 s=-1.D+00 call bwsm(nsim,vs2,c,t,s) sigma2=t esigma2=s c=1.D+00 t=-1.D+00 s=-1.D+00 c call bwsm(nsim,vsk,c,t,s) call SRTAD(vsk,1,nsim) eskew1=vsk(NINT(0.5*FLOAT(nsim)))-vsk(NINT(0.16*FLOAT(nsim))) eskew2=vsk(NINT(0.84*FLOAT(nsim)))-vsk(NINT(0.5*FLOAT(nsim))) skew=vsk(NINT(0.5*FLOAT(nsim))) c write(6,*) vsk(NINT(0.5*FLOAT(nsim))),eskew1,eskew2 c stop c eskew=s return end c c subroutine bwsm(n,xx,cc,theta,sigma) double precision par(10) double precision FUN,A,B,EPS,ANS,ERREST c double precision xx(*),x(n),cc ,theta,sigma,t1,s1,t2,s2 integer IROLD,IRNEW c common par external FUN,median,std if (n.eq.0) then theta=-1.D+00 sigma=-1.D+00 endif if (n.eq.1) then theta=x(1) sigma=0.D+00 endif kc=0 do k=1,n if (xx(k).eq.xx(1)) then kc=kc+1 endif enddo if (kc.eq.n) then theta=xx(1) sigma=0.1D+00*DAbs(theta) endif if ((n.le.1).or.(kc.eq.n)) then return endif do i=1,n x(i)=xx(i) enddo IRNEW=1 call ENTSRC(IROLD,IRNEW) par(1)=cc EPS=1.0D-10 A=-20.D+00 B=20.D+00 call DQUAD(FUN,A,B,EPS,ANS,ERREST) call ERROFF par(2)=ANS c write(6,*) ANS,ERREST call median(n,x,t1) call std(n,x,t1,s1) c write(6,*) 'inizio',t1,s1 k=0 101 k=k+1 s2=0.D+00 do i=1,n s2=s2+chi((x(i)-t1)/s1) enddo s2=DSqRt(s2)*s1/DSqRt(par(2)*DFLOAT(n-1)) t2=0.D+00 do i=1,n t2=t2+psi((x(i)-t1)/s2) enddo t2=t1+t2*s2/DFLOAT(n) if (((DAbs((t1-t2)/t1).gt.0.0001).or.(DAbs((s1-s2)/s1).gt.0.0001) $ ))then t1=t2 s1=s2 c write(6,*) k if (k.ge.200) goto 106 goto 101 endif 106 theta=t1 sigma=s1 103 return end c double precision function FUN(x) double precision par(10),x,chi common par FUN=chi(x)*1.0D+00/DSqRt(2.D+00*3.14159D+00)*DExp(-x**2/2.D+00) return end c double precision function psi(t) double precision par(10) double precision t common par c if (DAbs(t).le.1.0D+00) then psi=t*(1.D+00-t**2)**2 else psi=0.D+00 endif c c if ((t.gt.-1.D+00).and.(t.lt.0.D+00)) then c psi=t*(1.D+00-t**2)**2 c else c psi=0.D+00 c endif c if (t.ge.0.D+00) then c psi=t c endif return end c double precision function chi(t) double precision par(10) double precision t common par c if (DAbs(t).le.par(1)) then chi=t**2/2.D+00 else chi=par(1)**2/2.D+00 endif c if ((t.gt.-par(1))) then c chi=t**2/2.D+00 c else c chi=par(1)**2/2.D+00 c endif c return end c c subroutine stat(kk,vv,skew) double precision vv(*),skew,eskew,m,s m=0.D+00 do i=1,kk m=m+vv(i)/DFLOAT(kk) enddo s=0.D+00 do i=1,kk s=s+(vv(i)-m)**2 enddo s=dsqrt(s/DFLOAT(kk-1)) skew=0.D+00 do i=1,kk skew=skew+(vv(i)-m)**3/DFLOAT(kk) enddo skew=skew/s**3 return end c subroutine median(n,x,xmed) integer*4 n double precision x(*) integer*4 np, nm, nskip, j, niter, punt double precision xx, xp, xm, sumx, sum, eps, stemp, dum double precision ap, am, aa, a,xmed data big/1.0D30/, afac/1.5D+00/, amp/1.5D+00/, punt/9/ niter = 0 if(n .eq. 1) then xmed = x(1) return else if(n .eq. 2) then xmed = (x(1) + x(2))/2.0D+00 return endif a = 0.5D+00*(x(1)+x(n)) eps = abs(x(n)-x(1)) am = -big ap = big do while (.true.) niter = niter + 1 if(niter .gt. 50) then write(punt, '(a, i5, a)') 'median failed to converge for ', * n, ' points' write(punt, '(a, (e20.6))') 'data values were : ', * (x(i), i = 1, n) write(punt, '(a, e20.6)') 'leave median set at ', a xmed = a return endif sumx = 0.0 sum = 0.0 np = 0 nm = 0 nskip = 0 xm = -big xp = big j = 0 do while (j .lt. n) j = j + 1 xx = x(j) if (xx .ne. a) then if (xx .gt. a) then np = np + 1 if (xx .lt. xp) xp = xx else nm = nm + 1 if (xx .gt. xm) xm = xx endif dum = 1.0/(eps+abs(xx-a)) sum = sum + dum sumx = sumx + xx*dum else nskip = nskip + 1 endif enddo stemp = (sumx/sum) - a if (np-nm .ge. 2 .and. np-nm-nskip .ge. 2) then c guess too low am = a if ((stemp .le. 0.0) .or. (niter .gt. 25)) then aa = xp else aa = xp + stemp*amp endif if (aa .gt. ap) aa = 0.5*(a+ap) eps = afac*abs(aa-a) a = aa else if (nm-np .ge. 2 .and. nm-np-nskip .ge. 2) then c guess too high ap = a if ((stemp .gt. 0.0) .or. (niter .gt. 25)) then aa = xm else aa = xm + stemp*amp endif if (aa .lt. am) aa = 0.5*(a+am) eps = afac*abs(aa-a) a = aa else c success if ((n/2)*2 .eq. n) then c n is even if (nskip .gt. 1) then xmed = a else if (np .eq. nm) then xmed = 0.5*(xp+xm) else if (np .gt. nm) then xmed = 0.5*(a+xp) else xmed = 0.5*(xm+a) endif else c n is odd if (nskip .gt. 0) then xmed = a else if (np .eq. nm) then xmed = a else if (np .gt. nm) then xmed = xp else xmed = xm endif endif return endif enddo end c c subroutine medstd2(nn,xv,t1,s1) c double precision xv(*),xvv(10000000),xv2(10000000),x(10000000),t,s $ ,t1,s1,tt,ss,CCC,et1,et2 integer k1,k2,ko,ifail external statm,bws,median,median2 frac=2. do i=1,nn xvv(i)=xv(i) enddo k=0 300 k=k+1 if (k.eq.1) then call statm(nn,xvv,tt,ss) call median(nn,xvv,t) k2=nn endif if (k.gt.1) then t1=tt s1=ss k1=k2 k2=0 do j=1,nn if (DAbs(xvv(j)-tt).le.frac*ss) then k2=k2+1 xv2(k2)=xvv(j) endif enddo call statm(k2,xv2,tt,ss) call median(k2,xv2,tt) c write(6,FMT='(2I5,4F10.4)') k,k2,t1,s1,tt,ss if (((DAbs(ss-s1)/s1.lt.0.0457)).or.(k2.lt.nn/2)) goto310 endif goto 300 310 t=t1 k2=0 do j=1,nn if (DAbs(xvv(j)-tt).le.frac*s1) then k2=k2+1 xv2(k2)=xvv(j) endif enddo c call median2(k2,xv2,tt,et1,et2) c t=tt c call median(nn,xvv,t) c CCC=2.0D+00 c ifail=0 c t=-1.D+00 c s1=-1.D+00 c call bws(nn,xvv,t,s1,CCC,ifail) c s=1.25*s1/SqRt(FLOAT(k1)) c write(6,*) nn,k2 c return end c c subroutine statm(n,x,t,s) double precision x(10000000),t,s integer n t=0.D+00 do i=1,n t=t+x(i) enddo t=t/DFLOAT(n) s=0.D+00 do i=1,n s=s+(x(i)-t)**2 enddo s=SqRt(s/FLOAT(n-1)) return end c subroutine std(n,x,t1,s1) double precision x(*),t1,s1 c s1=0.D+00 do i=1,n s1=s1+DAbs(x(i)-t1)/DFLOAT(n)/0.798561862D+00 enddo c return end c