program addgal4 * compile using make parameter (itamax=15) parameter (nadd=40) parameter (pi=3.141593) include "setpar" real image(NDIM1,NDIM2) real scale /0.205/ common /imblk/ image integer x(nadd),y(nadd),id(nadd) integer ix,iy real mag ,magc real msurf ,msurfc real mpeak ,mpeakc real magdes,msurfdes real mtlist(100),mslist(100) real all(-15:15,-15:15,4,100) common /stack/ mtlist,mslist,all character*20 gimage(4) /'b.fits', & 'v.fits', & 'r.fits', & 'i.fits'/ real galr(-itamax:itamax,-itamax:itamax) real galq(-itamax:itamax,-itamax:itamax) real gal(-itamax:itamax,-itamax:itamax) real galcon(-itamax:itamax,-itamax:itamax) do 1 i=1,4 1 zp(i)=zp(i)+2.5*log10(et(i)) ! correct zero-points for exposure time magdes=23. msurfdes=23. idum=10 ! these are default values of course seeing= .9 / scale call condargr(1,magdes) call condargr(2,msurfdes) call condargi(3,idum) print*,'loading real galaxies' call loadgals() print*,'loading real galaxies: DONE' * print*,'generating galaxy' * call gengal(ceni,tau,seeing,gal,galcon) * print*,'generating galaxy: DONE' * generate galaxy locations and galaxy id's do 5 i=1,nadd x(i)=nint( (ran3(idum)*.9+.05) * real(ndim1) ) y(i)=nint( (ran3(idum)*.9+.05) * real(ndim2) ) 4 id(i)=int( ran3(idum)*100.)+1 if (mtlist(id(i)).gt.magdes) goto 4 5 continue * big loop over all colours do 100 ic=1,4 call movein(orient//gimage(ic),image,NDIM1,NDIM2) count=0. do 10 i=1,nadd ix=x(i) iy=y(i) * ix=mod(100*(i-1),900)+100 * iy=int((i-1)/9)*100+100 * id(i)=i count=count+1. call pickgal(galr,ic,id(i)) call tovals(mtlist(id(i)),mslist(id(i)),magdes,msurfdes,side,up) call smear(galr,side,up,galq) call insgal(ix,iy,galq) if (ic.eq.4) write(21,*) ix,iy,mtlist(id(i)),mslist(id(i)),id(i) 10 continue print*,'moving out td'//gimage(ic) call writefitsreal('td'//gimage(ic),image,NDIM1,NDIM2) 100 continue ic=4 write(*,*) ' mag mpeak' write(*,'(2f8.3)') magdes,msurfdes 9998 format (2i5,3f8.3,2f10.5) 9999 format(2f8.2,4(f9.4,f7.4)) stop end *********************** *********************** subroutine tovals(magin,msurfin,magout,msurfout,side,up) real magin,msurfin,magout,msurfout up=10**(0.4*(magin-magout)) side=10**(0.2*(magin-magout-msurfin+msurfout)) return end *********************** subroutine smear(galin,side,up,galout) parameter (itamax=15) real galin(-itamax:itamax,-itamax:itamax) real galout(-itamax:itamax,-itamax:itamax) do 10 j=-itamax,itamax do 10 i=-itamax,itamax x=real(i)/side y=real(j)/side if (max(abs(x),abs(y)).gt.itamax) then val=0. goto 9 endif i1=ibot(x) i2=itop(x) j1=ibot(y) j2=itop(y) t=x-i1 u=y-j1 val= (1.-t)*(1.-u) * galin(i1,j1) val=val+(1.-t)*( u) * galin(i1,j2) val=val+( t)*(1.-u) * galin(i2,j1) val=val+( t)*( u) * galin(i2,j2) 9 continue val=val/side/side*up 10 galout(i,j)=val return end *********************** subroutine getmag(gal,zp,scale,mag,msurf) parameter (itamax=15) real gal(-itamax:itamax,-itamax:itamax) real mag,msurf flux=0. fluxmax=gal(0,0) do 10 j=-itamax,itamax do 10 i=-itamax,itamax flux=flux+gal(i,j) 10 fluxmax=max(fluxmax,gal(i,j)) fluxmax=fluxmax/scale/scale mag = -2.5*log10(flux )+zp msurf= -2.5*log10(fluxmax)+zp return end *********************** subroutine pickgal(galr,ic,id) parameter (itamax=15) real mtlist(100),mslist(100) real all(-itamax:itamax,-itamax:itamax,4,100) common /stack/ mtlist,mslist,all real galr(-itamax:itamax,-itamax:itamax) do 10 j=-itamax,itamax do 10 i=-itamax,itamax 10 galr(i,j)=all(i,j,ic,id) return end *********************** subroutine loadgals() parameter (itamax=15) real mtlist(100),mslist(100) real all(-itamax:itamax,-itamax:itamax,4,100) common /stack/ mtlist,mslist,all open(31,file='mmstack',form='unformatted') read(31) mtlist read(31) mslist read(31) all close(31) return end *********************** subroutine putgal(ixc,iyc,cen,tau,flux,surf,pix,pixq,sigsky) PARAMETER (NDIM1=1024, NDIM2=1024) real image(NDIM1,NDIM2) common /imblk/ image ita=nint(tau*5.) ita=max(ita,15) clip=0.70000*sigsky ! a very low threshold if you ask me flux=0. surf=0. pix=0. do 10 ix=-ita,ita do 10 iy=-ita,ita ixq=ixc+ix iyq=iyc+iy r=sqrt(real(ix*ix+iy*iy)) fluxq=cen*exp(-r/tau) flux=flux+fluxq if (fluxq.gt.clip) then pix=pix+1. surf=surf+fluxq endif if (ixq.lt.1) goto 9 if (iyq.lt.1) goto 9 if (ixq.gt.ndim1) goto 9 if (iyq.gt.ndim2) goto 9 image(ixq,iyq)=image(ixq,iyq)+fluxq 9 continue 10 continue rq=-log(7.78415/cen)*tau pixq=3.14159*rq*rq surf=surf/pix return end *********************** subroutine gengal(cen,tau,seeing,gal,galcon) parameter (itamax=15) parameter (pi=3.141593) real gal(-itamax:itamax,-itamax:itamax) real galcon(-itamax:itamax,-itamax:itamax) ita=nint(tau*5.) ita=min(ita,itamax) print*,' generating disk' * generate the exponential disk do 10 ix=-ita,ita do 10 iy=-ita,ita r=sqrt(real(ix*ix+iy*iy)) 10 gal(ix,iy)=cen*exp(-r/tau) print*,' done' * convolve with seeing pixel radius Gaussian * seeing is the FWHM print*,' doing the seeing' sigma=0.6006*seeing ! see DET.log 22/10/2001 for FHWM -->sigma irel = int (4 * sigma) ! the radius out which we should integrate print*,' sigma (pixels):', sigma if (irel.lt.1) PAUSE 'gengal: write more code' anorm = 1. / (pi* sigma*sigma ) ! see DET.log 22/10/2001 do 20 ix=-ita,ita do 20 iy=-ita,ita fluxq=0. do 30 jx=-irel,irel do 30 jy=-irel,irel kx=jx+ix ky=jy+iy if (max(abs(kx),abs(ky)).le.ita) then r=sqrt( real(jx*jx + jy*jy) ) gaussq= anorm * exp ( - ((r/sigma)**2) ) fluxq=fluxq+gaussq*gal(kx,ky) endif 30 continue 20 galcon(ix,iy)=fluxq print*,' done' return end *********************** subroutine insgal(ixc,iyc,gal) parameter (itamax=15) PARAMETER (NDIM1=1024, NDIM2=1024) real gal(-itamax:itamax,-itamax:itamax) real image(NDIM1,NDIM2) common /imblk/ image ita=max(ita,itamax) do 10 ix=-itamax,itamax do 10 iy=-itamax,itamax ixq=ixc+ix iyq=iyc+iy if (ixq.lt.1) goto 9 if (iyq.lt.1) goto 9 if (ixq.gt.ndim1) goto 9 if (iyq.gt.ndim2) goto 9 image(ixq,iyq)=image(ixq,iyq)+gal(ix,iy) 9 continue 10 continue return end ***********************