program addgal
*  compile using make     
	PARAMETER (NDIM1=1024, NDIM2=1024)
	parameter (itamax=15)
	parameter (pi=3.141593)
      real image(NDIM1,NDIM2)
	common /imblk/ image
      real zp /29.81/
      real et /1./
      real sigsky /1.63675/
      real scale /0.205/
      real x,y
	real mag    ,magc
	real msurf  ,msurfc
	real mpeak  ,mpeakc
        real gal(-itamax:itamax,-itamax:itamax)
	real galcon(-itamax:itamax,-itamax:itamax)

	zp=zp+2.5*log10(et)
	
	call movein('ti.fits',image,NDIM1,NDIM2)
	cen=50. ! mag scales (and msurf, basically) as -2.5*log10(cen)
	tau=5.  ! msurf is independent of tau, mag scales as -5*log10(tau)
	idum=10
	seeing= .9 / scale
		
	call condargr(1,cen)
	call condargr(2,tau)
	call condargi(3,idum)

	print*,'generating galaxy'
	call gengal(cen,tau,seeing,gal,galcon)
	print*,'generating galaxy: DONE'
		
	do 10 i=1,40
	ix=nint(ran3(idum)*real(ndim1)*.9+.05)  !!!! there is bug here
	iy=nint(ran3(idum)*real(ndim2)*.9+.05)
	count=count+1.
*	call putgal(ix,iy,cen,tau,flux,surf,pix,pixq,sigsky)
	call insgal(ix,iy,galcon)

	fluxc=2.*pi*cen*tau*tau
	mag=-2.5*log10(fluxc)+zp
	mpeak=-2.5*log10(cen/scale/scale)+zp
	musrf=0.   ! OK actually it isn't
	
9998	format (2i5,3f8.3,2f10.5)
	write(21,9998) ix,iy,mag,msurf,mpeak,cen,tau
10	continue

	print*,'moving out'
	call writefitsreal('tdone.fits',image,NDIM1,NDIM2)
	
	print*,mag,mpeak
	write(*,'(2f8.3)') mag,mpeak

9999	format(2f8.2,4(f9.4,f7.4))
	stop
      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
***********************	
	
	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 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
***********************