diff -r -u -b -P SFH-1.1/grid/mkgrid.f SFH-1.1-p1/grid/mkgrid.f --- SFH-1.1/grid/mkgrid.f 2004-05-13 12:59:57.000000000 -0700 +++ SFH-1.1-p1/grid/mkgrid.f 2004-06-02 11:09:04.000000000 -0700 @@ -8,7 +8,7 @@ nbin = 5 do icmd=1,ncmd - open(unit=11,file="grid"//icmd+48//".index",type="unknown") + open(unit=11,file="grid"//char(icmd+48)//".index",type="unknown") do ix=1,nx(icmd) do iy=1,ny Only in SFH-1.1/input/crowd: test.cfn diff -r -u -b -P SFH-1.1/repop/Makefile SFH-1.1-p1/repop/Makefile --- SFH-1.1/repop/Makefile 2004-05-07 12:20:44.000000000 -0700 +++ SFH-1.1-p1/repop/Makefile 2004-06-03 11:17:08.000000000 -0700 @@ -4,8 +4,9 @@ SYNTHDIR = ../synthcode COMDIR = ../commoncode OBJECTS = repop.o axes.o bin_dmags.o bin_interp.o bincombine.o \ -crowd.o cumfcn.o dmags.o fakecrowd.o mkcfn.o outpix.o pickphot.o \ -readcrowd.o readphot.o redden.o random.o parseline.o utilities.o +crowd.o cumfcn.o delta_axes.o dmags.o fakecrowd.o mkcfn.o outpix.o \ +pickphot.o readcrowd.o readphot.o redden.o random.o parseline.o \ +utilities.o repop : $(OBJECTS) $(SYNTHDIR)/synth.h $(FORT) -o repop $(OBJECTS) @@ -31,6 +32,9 @@ cumfcn.o : $(SYNTHDIR)/cumfcn.f $(FORT) $(CFLAGS) $(SYNTHDIR)/cumfcn.f +delta_axes.o : $(SYNTHDIR)/delta_axes.f + $(FORT) $(CFLAGS) $(SYNTHDIR)/delta_axes.f + dmags.o : $(SYNTHDIR)/dmags.f $(FORT) $(CFLAGS) $(SYNTHDIR)/dmags.f diff -r -u -b -P SFH-1.1/repop/repop.f SFH-1.1-p1/repop/repop.f --- SFH-1.1/repop/repop.f 2004-05-19 17:30:40.000000000 -0700 +++ SFH-1.1-p1/repop/repop.f 2004-06-03 11:17:22.000000000 -0700 @@ -22,6 +22,7 @@ integer lockflag,i,ib,oldb,idrop(MCMD) integer jnext,jpercent,iseed integer nstars,np + integer err_method, interp_errs real age,t,z,t2,z2,dtime real x(MCMD),y(MCMD),x0(MCMD),y0(MCMD) real xmax(MCMD),ycmax(MCMD),ypmax(MCMD) @@ -50,15 +51,12 @@ character*40 crowd1,crowd2 common /c_cfile/ crowd1,crowd2 - integer nx(MCMD),ny(MCMD),interp_errs, fake_errs - real dbinx,dbiny,dbiny2 - common /c_crowd/ interp_errs,fake_errs,dbinx,dbiny,nx,ny - real emin,emax,epix common /c_chist/ emin,emax,epix - real xmin(MCMD),ymin(MCMD) - common /c_min/ xmin,ymin + integer nx(MCMD),ny(MCMD) + real dbinx,dbiny,xmin(MCMD),ymin(MCMD) + common /c_crowd/ dbinx,dbiny,nx,ny,xmin,ymin character*8 xeqn(MCMD),yeqn(MCMD) common /c_eqn/ xeqn,yeqn @@ -118,7 +116,7 @@ c *** Miscellaneous *** call fetchint( 1, iverb ) call fetchint( 1, interp_errs ) - call fetchint( 1, fake_errs ) + call fetchint( 1, err_method ) call fetchint( 1, nscale ) call fetchint( 1, iseed ) call fetchreal( 1, mass1 ) @@ -271,7 +269,7 @@ end do ! Add photometric scatter - if ( fake_crowd .eq. 1 ) then + if ( err_method .eq. 1 ) then call fakecrowd( x, y, idrop ) else call crowd( x, y, idrop ) @@ -302,8 +300,6 @@ close(20+icmd) end do - 4 format(f5.2,2(3x,a24)) - 5 format(i5,f6.2,2(3x,a24)) 7 format(4f7.3) stop diff -r -u -b -P SFH-1.1/sfhcode/dchi.f SFH-1.1-p1/sfhcode/dchi.f --- SFH-1.1/sfhcode/dchi.f 2004-05-18 17:54:16.000000000 -0700 +++ SFH-1.1-p1/sfhcode/dchi.f 2004-06-02 11:23:06.000000000 -0700 @@ -47,7 +47,7 @@ do 11 j=1,JMAX dx = dx*0.5 xmid = dchi + dx - fmid = gamma_q( hdof, dble(xmid)/2.0D0 ) + fmid = real(gamma_q( hdof, dble(xmid)/2.0D0 )) - q if (fmid.le.0.) dchi = xmid if (abs(dx).lt.xacc .or. fmid.eq.0.) return 11 continue Only in SFH-1.1/sfhcode: test_dchi Only in SFH-1.1/sfhcode: test_dchi.f Only in SFH-1.1/sfhcode: test_dchi.o diff -r -u -b -P SFH-1.1/synthcode/bincombine.f SFH-1.1-p1/synthcode/bincombine.f --- SFH-1.1/synthcode/bincombine.f 2004-05-11 11:32:14.000000000 -0700 +++ SFH-1.1-p1/synthcode/bincombine.f 2004-06-02 13:59:55.000000000 -0700 @@ -29,9 +29,9 @@ integer ndrop(BINX,BINY),nstars(BINX,BINY) real dx(2,BINX,BINY,MAXPOP),dx2(2,BINX,BINY,MAXPOP) - integer nx(MCMD),ny(MCMD),interp_errs,fake_errs - real dbinx,dbiny - common /c_crowd/ interp_errs,fake_errs,dbinx,dbiny,nx,ny + integer nx(MCMD),ny(MCMD) + real dbinx,dbiny,xmin(MCMD),ymin(MCMD) + common /c_crowd/ dbinx,dbiny,nx,ny,xmin,ymin c Copy dx to dx2 do i=1,2 diff -r -u -b -P SFH-1.1/synthcode/crowd.f SFH-1.1-p1/synthcode/crowd.f --- SFH-1.1/synthcode/crowd.f 2004-05-11 11:32:59.000000000 -0700 +++ SFH-1.1-p1/synthcode/crowd.f 2004-06-02 14:07:43.000000000 -0700 @@ -128,12 +128,9 @@ integer nmag,ncmd common /c_nmag/ nmag,ncmd - integer nx(MCMD),ny(MCMD), interp_errs, fake_errs - real dbinx,dbiny - common /c_crowd/ interp_errs,fake_errs,dbinx,dbiny,nx,ny - - real xmin(MCMD),ymin(MCMD) - common /c_min/ xmin,ymin + integer nx(MCMD),ny(MCMD) + real dbinx,dbiny,xmin(MCMD),ymin(MCMD) + common /c_crowd/ dbinx,dbiny,nx,ny,xmin,ymin real pdrop(MCMD,BINX,BINY) real cumerr(MCMD,BINX,BINY,NBINS*NBINS) @@ -169,7 +166,7 @@ goto 200 ! skip phot errs stuff endif -ccc Uncomment to interpolate crowding errors between CMD bins + if ( interp_errs.eq.1 ) then ! pick bracketing crowding bins for bin-interpolation jx = ix + 1 jy = iy + 1 @@ -203,22 +200,20 @@ ! The interpolated values are returned through (dx,dy) and p0. call bin_interp( icmd, ix, iy, dxbin, dybin, x dm, pdrop, dx, dy, p0 ) -ccc (end of block to uncomment for bin-interpolation) + else + ! Choose a random number for indexing the cumulative + ! error functions. Because we are recasting double precision + ! to real, it is possible to have a value that is equal + ! to 1.0 (the original double was less than 1.0). When + ! this occurs, draw a different random number. + 199 call random( dr, 0.0D0, 1.0D0 ) + rr = real(dr) + if ( rr.eq.1.0 ) goto 199 + + p0 = pdrop(icmd,ix,iy) ! just use the bin's dropout fraction -ccc Uncomment this block if NOT interpolating between CMD bins -! ! Choose a random number for indexing the cumulative -! ! error functions. Because we are recasting double precision -! ! to real, it is possible to have a value that is equal -! ! to 1.0 (the original double was less than 1.0). When -! ! this occurs, draw a different random number. -! 199 call random( dr, 0.0D0, 1.0D0 ) -! rr = real(dr) -! if ( rr.eq.1.0 ) goto 199 -! -! p0 = pdrop(icmd,ix,iy) ! just use the bin's dropout fraction -! -! call dmags( icmd, ix, iy, rr, dx, dy ) -ccc (end of block to be commented-out for bin interpolation) + call dmags( icmd, ix, iy, rr, dx, dy ) + endif ! Now that we have a dropout fraction, we can test for dropout call random(f,0.0D0,1.0D0) diff -r -u -b -P SFH-1.1/synthcode/delta_axes.f SFH-1.1-p1/synthcode/delta_axes.f --- SFH-1.1/synthcode/delta_axes.f 1969-12-31 17:00:00.000000000 -0700 +++ SFH-1.1-p1/synthcode/delta_axes.f 2004-06-02 13:46:41.000000000 -0700 @@ -0,0 +1,104 @@ + subroutine delta_axes(dmag,dx,dy) + +c This file is part of StarFISH +c (C) 2004 by Jason Harris +c +c -=> StarFISH is free software; you can redistribute it +c -=> and/or modify it under the terms of the GNU General Public +c -=> License (GPL) as published by the Free Software Foundation; +c -=> either version 2 of the License, or (at your option) any +c -=> later version. +c +c Constructs CMD offsets according to the equation strings +c (xeqn(ncmd), yeqn(ncmd)) found in synth.dat. +c +c This is just like axes(), but it is used to determine the +c delta-mags in the X- and Y-directions. This is a special +c case because the delta-mags also encode whether the star +c is a dropout. + + include 'synth.h' + + character*8 str + integer icmd,xy,col(MCMD,2,2),iflag + real dx(MCMD),dy(MCMD),sum + real dmag(MMAG) + logical number,add(MCMD,2) +c number = .true. : expect a number +c number = .false. : expect +/- + + save add, col, iflag + + character*8 xeqn(MCMD),yeqn(MCMD) + common /c_eqn/ xeqn, yeqn + + integer nmag,ncmd + common /c_nmag/ nmag,ncmd + +c On first call, parse the eqn strings and store numbers, operators +c in ixcol/iycol, xadd/yadd + if (iflag.ne.1) then + iflag = 1 + + do icmd=1,ncmd + do xy=1,2 + if (xy.eq.1) str = xeqn(icmd) + if (xy.eq.2) str = yeqn(icmd) + jcol = 1 + col(icmd,xy,2) = 0 ! default to no second argument + number = .true. + + do i=1,len(str) + if (str(i:i).ne.' ') then + value = ichar(str(i:i)) - 48 + if (number) then + if (value.lt.0.or.value.gt.9) then ! not a digit + write(*,*) str + pause 'parse error: digit expected' + else + col(icmd,xy,jcol) = value + jcol = jcol + 1 + if (jcol.eq.3) goto 10 ! done + endif + else + add(icmd,xy) = .true. + if (value.eq.-3) add(icmd,xy) = .false. ! read '-' + endif + number = .not. number + endif + end do + 10 continue + end do + end do + endif + + do icmd=1,ncmd + do xy=1,2 + if (col(icmd,xy,1).eq.0) + x pause 'Error parsing magnitudes in axes.f ' + + ! If either mag is equal to 99.9, the star is a dropout. + ! Instead of summing/subtracting the two mags in this + ! case, leave dx/dy as 99.9 so the fact that it dropped + ! is preserved. + if ( dmag(col(icmd,xy,1)).gt.90.0 .or. + x dmag(col(icmd,xy,2)).gt.90.0 ) then + sum = 99.9 + else + sum = dmag(col(icmd,xy,1)) + if (col(icmd,xy,2).ne.0) then + if (add(icmd,xy)) then + sum = sum + dmag(col(icmd,xy,2)) + else + sum = sum - dmag(col(icmd,xy,2)) + endif + endif + endif + + if (xy.eq.1) dx(icmd) = sum + if (xy.eq.2) dy(icmd) = sum + end do + end do + + return + end diff -r -u -b -P SFH-1.1/synthcode/dmags.f SFH-1.1-p1/synthcode/dmags.f --- SFH-1.1/synthcode/dmags.f 2004-05-11 11:33:37.000000000 -0700 +++ SFH-1.1-p1/synthcode/dmags.f 2004-06-02 11:14:27.000000000 -0700 @@ -40,13 +40,19 @@ endif ! find the index kk at which the cumfcn crosses the number rr - kk = 0 - do k=1,nbin*nbin - if ( cumerr(icmd,ix,iy,k).gt.rr ) then - kk = k - goto 10 + dk = int(nbin*nbin/2) + kk = dk + + 10 c = cumerr(icmd,ix,iy,kk) + + dk = int(dk/2) + if ( dk.le.0 ) goto 20 + if ( c.gt.rr ) then + kk = kk - dk + else + kk = kk + dk endif - end do + goto 10 write(*,*) 'Error: point not found in cumulative err function!' write(*,*) icmd,ix,iy,rr,cumerr(icmd,ix,iy,nbin*nbin) @@ -57,7 +63,7 @@ ! Convert from the 1-D cum.fcn. index (kk) to ! the 2-D delta-mag pixel indices (kx,ky) - 10 ky = int( kk/nbin ) + 1 + 20 ky = int( kk/nbin ) + 1 kx = kk - nbin*( ky - 1 ) ! determine the values of the bracketing cumfcn bins (x1 and x2) @@ -92,7 +98,7 @@ dx = xb( kx-1 ) + f*( xb( kx ) - xb( kx-1 ) ) ! the dY value is the row's bin value, plus random scatter. call random( dr, -0.50D0, 0.50D0 ) - dy = real(dr) + xb( ky ) + dy = epix*real(dr) + xb( ky ) return end diff -r -u -b -P SFH-1.1/synthcode/Makefile SFH-1.1-p1/synthcode/Makefile --- SFH-1.1/synthcode/Makefile 2004-05-07 11:55:49.000000000 -0700 +++ SFH-1.1-p1/synthcode/Makefile 2004-06-03 11:15:18.000000000 -0700 @@ -4,8 +4,9 @@ COMDIR = ../commoncode BINDIR = ../. OBJECTS = synth.o axes.o bin_dmags.o bin_interp.o bincombine.o \ -crowd.o cumfcn.o dmags.o fakecrowd.o mkcfn.o outpix.o pickphot.o \ -readcrowd.o readphot.o redden.o random.o parseline.o utilities.o +crowd.o cumfcn.o delta_axes.o dmags.o fakecrowd.o mkcfn.o outpix.o \ +pickphot.o readcrowd.o readphot.o redden.o random.o parseline.o \ +utilities.o synth : $(OBJECTS) synth.h $(FORT) -o synth $(OBJECTS) diff -r -u -b -P SFH-1.1/synthcode/mkcfn.f SFH-1.1-p1/synthcode/mkcfn.f --- SFH-1.1/synthcode/mkcfn.f 2004-05-12 17:34:54.000000000 -0700 +++ SFH-1.1-p1/synthcode/mkcfn.f 2004-06-02 13:58:59.000000000 -0700 @@ -32,12 +32,9 @@ character*40 crowd1,crowd2 common /c_cfile/ crowd1,crowd2 - integer nx(MCMD),ny(MCMD),interp_errs,fake_errs - real dbinx,dbiny - common /c_crowd/ interp_errs,fake_errs,dbinx,dbiny,nx,ny - - real xmin(MCMD),ymin(MCMD) - common /c_min/ xmin,ymin + integer nx(MCMD),ny(MCMD) + real dbinx,dbiny,xmin(MCMD),ymin(MCMD) + common /c_crowd/ dbinx,dbiny,nx,ny,xmin,ymin real emin,emax,dbin1 common /c_chist/ emin,emax,dbin1 @@ -89,7 +86,7 @@ c convert mag and dmag into color/mag CMD positions (x/y) call axes(mag,x,y) - call axes(dmag,xoff,yoff) + call delta_axes(dmag,xoff,yoff) c Determine to which crowding bin the current star belongs. ix = int((x(icmd) - xmin(icmd))/dbinx) + 1 diff -r -u -b -P SFH-1.1/synthcode/readcrowd.f SFH-1.1-p1/synthcode/readcrowd.f --- SFH-1.1/synthcode/readcrowd.f 2004-05-12 17:57:58.000000000 -0700 +++ SFH-1.1-p1/synthcode/readcrowd.f 2004-06-02 13:59:28.000000000 -0700 @@ -22,9 +22,9 @@ character*40 crowd1,crowd2 common /c_cfile/ crowd1,crowd2 - integer nx(MCMD),ny(MCMD),interp_errs,fake_errs - real dbinx,dbiny - common /c_crowd/ interp_errs,fake_errs,dbinx,dbiny,nx,ny + integer nx(MCMD),ny(MCMD) + real dbinx,dbiny,xmin(MCMD),ymin(MCMD) + common /c_crowd/ dbinx,dbiny,nx,ny,xmin,ymin real emin,emax,dbin1 common /c_chist/ emin,emax,dbin1 diff -r -u -b -P SFH-1.1/synthcode/readphot.f SFH-1.1-p1/synthcode/readphot.f --- SFH-1.1/synthcode/readphot.f 2004-05-19 15:06:56.000000000 -0700 +++ SFH-1.1-p1/synthcode/readphot.f 2004-06-02 14:28:47.000000000 -0700 @@ -14,7 +14,7 @@ c include 'synth.h' - character*24 libfile + character*40 libfile integer n,j real f diff -r -u -b -P SFH-1.1/synthcode/redden.f SFH-1.1-p1/synthcode/redden.f --- SFH-1.1/synthcode/redden.f 2004-05-11 11:43:55.000000000 -0700 +++ SFH-1.1-p1/synthcode/redden.f 2004-06-02 14:28:55.000000000 -0700 @@ -12,7 +12,7 @@ c include 'synth.h' - character*24 avfile + character*40 avfile integer i,n,ntot real av(NAV) diff -r -u -b -P SFH-1.1/synthcode/synth.f SFH-1.1-p1/synthcode/synth.f --- SFH-1.1/synthcode/synth.f 2004-05-13 11:38:02.000000000 -0700 +++ SFH-1.1-p1/synthcode/synth.f 2004-06-03 11:15:58.000000000 -0700 @@ -30,6 +30,7 @@ integer i,redflag,ib,oldb,idrop(MCMD) integer jnext,jpercent integer nstars,nadd,np + integer interp_errs, err_method real x(MCMD),y(MCMD),age,dpix real xmax(MCMD),ycmax(MCMD),ypmax(MCMD) real gamma,mass1,mass2,mrange,fsum,frac @@ -57,16 +58,13 @@ character*40 crowd1,crowd2 common /c_cfile/ crowd1,crowd2 - integer nx(MCMD),ny(MCMD),interp_errs, fake_errs - real dbinx,dbiny,dbiny2 - common /c_crowd/ interp_errs,fake_errs,dbinx,dbiny,nx,ny + integer nx(MCMD),ny(MCMD) + real dbinx,dbiny,xmin(MCMD),ymin(MCMD) + common /c_crowd/ dbinx,dbiny,nx,ny,xmin,ymin real emin,emax,epix common /c_chist/ emin,emax,epix - real xmin(MCMD),ymin(MCMD) - common /c_min/ xmin,ymin - character*8 xeqn(MCMD),yeqn(MCMD) common /c_eqn/ xeqn,yeqn @@ -117,7 +115,7 @@ c *** Miscellaneous *** call fetchint( 5, iverb ) call fetchint( 5, interp_errs ) - call fetchint( 5, fake_errs ) + call fetchint( 5, err_method ) call fetchint( 5, nscale ) call fetchint( 5, iseed ) call fetchreal( 5, mass1 ) @@ -233,7 +231,7 @@ call axes(starmag,x,y) ! Add photometric scatter - if ( fake_errs .eq. 1 ) then + if ( err_method .eq. 1 ) then call fakecrowd( x, y, idrop ) else call crowd( x, y, idrop ) @@ -263,8 +261,5 @@ close(20+icmd) end do - 4 format(f5.2,2(3x,a24)) - 5 format(i5,f6.2,2(3x,a24)) - stop end diff -r -u -b -P SFH-1.1/testpop/Makefile SFH-1.1-p1/testpop/Makefile --- SFH-1.1/testpop/Makefile 2004-05-07 12:08:20.000000000 -0700 +++ SFH-1.1-p1/testpop/Makefile 2004-06-03 11:17:40.000000000 -0700 @@ -4,8 +4,9 @@ SYNTHDIR = ../synthcode COMDIR = ../commoncode OBJECTS = testpop.o axes.o bin_dmags.o bin_interp.o bincombine.o \ -crowd.o cumfcn.o dmags.o fakecrowd.o mkcfn.o outpix.o pickphot.o \ -readcrowd.o readphot.o redden.o random.o parseline.o utilities.o +crowd.o cumfcn.o delta_axes.o dmags.o fakecrowd.o mkcfn.o outpix.o \ +pickphot.o readcrowd.o readphot.o redden.o random.o parseline.o \ +utilities.o testpop : $(OBJECTS) $(SYNTHDIR)/synth.h $(FORT) -o testpop $(OBJECTS) @@ -31,6 +32,9 @@ cumfcn.o : $(SYNTHDIR)/cumfcn.f $(FORT) $(CFLAGS) $(SYNTHDIR)/cumfcn.f +delta_axes.o : $(SYNTHDIR)/delta_axes.f + $(FORT) $(CFLAGS) $(SYNTHDIR)/delta_axes.f + dmags.o : $(SYNTHDIR)/dmags.f $(FORT) $(CFLAGS) $(SYNTHDIR)/dmags.f diff -r -u -b -P SFH-1.1/testpop/testpop.f SFH-1.1-p1/testpop/testpop.f --- SFH-1.1/testpop/testpop.f 2004-05-19 15:10:32.000000000 -0700 +++ SFH-1.1-p1/testpop/testpop.f 2004-06-03 11:18:21.000000000 -0700 @@ -22,6 +22,7 @@ integer lockflag,sfrflag,i,ib,oldb,idrop(MCMD) integer jnext,jpercent,iseed integer nstars,np + integer err_method, interp_errs real age,t,z,t2,z2,dtime real x(MCMD),y(MCMD),x0(MCMD),y0(MCMD) real xmax(MCMD),ycmax(MCMD),ypmax(MCMD) @@ -50,16 +51,13 @@ character*40 crowd1,crowd2 common /c_cfile/ crowd1,crowd2 - integer nx(MCMD),ny(MCMD),interp_errs, fake_errs - real dbinx,dbiny,dbiny2 - common /c_crowd/ interp_errs,fake_errs,dbinx,dbiny,nx,ny + integer nx(MCMD),ny(MCMD) + real dbinx,dbiny,xmin(MCMD),ymin(MCMD) + common /c_crowd/ dbinx,dbiny,nx,ny,xmin,ymin real emin,emax,epix common /c_chist/ emin,emax,epix - real xmin(MCMD),ymin(MCMD) - common /c_min/ xmin,ymin - character*8 xeqn(MCMD),yeqn(MCMD) common /c_eqn/ xeqn,yeqn @@ -118,7 +116,7 @@ c *** Miscellaneous *** call fetchint( 1, iverb ) call fetchint( 1, interp_errs ) - call fetchint( 1, fake_errs ) + call fetchint( 1, err_method ) call fetchint( 1, nscale ) call fetchint( 1, iseed ) call fetchreal( 1, mass1 ) @@ -303,7 +301,7 @@ end do ! Add photometric scatter - if ( fake_crowd .eq. 1 ) then + if ( err_method .eq. 1 ) then call fakecrowd( x, y, idrop ) else call crowd( x, y, idrop )