
# set global parameters for the simulations
macro mkimages.init
 if ($0 != 7)
  echo "USAGE: mkimage.empty (Nx) (Ny) (RDNOISE) (SKY) (Nstars) (psfSigma)"
  break
 end

 $mkiNx = $1
 $mkiNy = $2
 $mkiRD = $3
 $mkiSKY = $4
 $mkiPSFsigma = $6

 gaussdev mki_noisev {$mkiNx*$mkiNy} 0.0 1.0
 dimenup mki_noisev mki_noisei $mkiNx $mkiNy
 mcreate blank $mkiNx $mkiNy
end

# generate a fake image with bias and sky, no stars
macro mkimage.empty
 if ($0 != 2)
  echo "USAGE: mkimage.empty (buffer)"
  break
 end

 set $1 = blank + $mkiRD*mki_noisei + sqrt($mkiSKY)*mki_noisei + $mkiSKY
end

# set global parameters for the simulations
macro mkimages.stars
 if ($0 != 7)
  echo "USAGE: mkimage.empty (Nx) (Ny) (RDNOISE) (SKY) (Nstars) (psfSigma)"
  break
 end

 $mkiNx = $1
 $mkiNy = $2
 $mkiRD = $3
 $mkiSKY = $4
 $mkiPSFsigma = $6

 gaussdev mki_noisev {$mkiNx*$mkiNy} 0.0 1.0
 dimenup mki_noisev mki_noisei $mkiNx $mkiNy
 mcreate blank $mkiNx $mkiNy

 # mkstars $5 0.6
 mkstars $5 0.4
 # mkstars $5 0.3
end

# generate a fake image with bias, sky, stars, gaussian PFS
macro mkimage.stars.gauss
 if ($0 != 2)
  echo "USAGE: mkimage.stars (buffer)"
  break
 end

 set _tmp.stars = mki_stars_image

 mkimage.smooth _tmp.stars _tmp.smooth $mkiPSFsigma mkgauss

 # noise sigma = sqrt (_tmp.smooth + $mkiRD^2 + $mkiSKY)
 set $1 = _tmp.smooth + $mkiSKY + sqrt(_tmp.smooth + $mkiSKY + $mkiRD^2)*mki_noisei
end

# generate a fake image with bias, sky, stars, gaussian PFS
macro mkimage.stars.power
 if ($0 != 2)
  echo "USAGE: mkimage.stars (buffer)"
  break
 end

 set _tmp.stars = mki_stars_image

 mkimage.smooth _tmp.stars _tmp.smooth $mkiPSFsigma mkpower

 # noise sigma = sqrt (_tmp.smooth + $mkiRD^2 + $mkiSKY)
 set $1 = _tmp.smooth + $mkiSKY + sqrt(_tmp.smooth + $mkiSKY + $mkiRD^2)*mki_noisei
end

# generate a fake image with bias, sky, stars, gaussian PFS
macro mkimage.stars.wing
 if ($0 != 2)
  echo "USAGE: mkimage.stars (buffer)"
  break
 end

 set _tmp.stars = mki_stars_image

 mkimage.smooth _tmp.stars _tmp.smooth $mkiPSFsigma mkwing

 # noise sigma = sqrt (_tmp.smooth + $mkiRD^2 + $mkiSKY)
 set $1 = _tmp.smooth + $mkiSKY + sqrt(_tmp.smooth + $mkiSKY + $mkiRD^2)*mki_noisei
end

# generate the vectors of star coordinates, fluxes and magnitudes
macro mkstars
 if ($0 != 3)
   echo "USAGE: mkstars (Nstars) (alpha)"
   break
 end

 local sigma alpha Nstars
 $Nstars = $1
 $alpha = $2
 # alpha is dlogN/dmag = 0.6 for N ~ S^-3/2

 # $sigma is flux of star with S/N = 1 (sigma_sky * 4pi sigma_psf^2)
 $sigma = sqrt($mkiRD^2 + $mkiSKY) * 4*3.1416*$mkiPSFsigma^2
 
 # faintest magnitudes:
 # include stars down to 1 sigma
 $M1 = -2.5*log($sigma)

 # brightest magnitude in sim
 $M0 = $M1 - log($alpha*$Nstars/0.43)/$alpha

 create _tmp 0 $Nstars
 set mki_stars_x = $mkiNx * rnd(_tmp)
 set mki_stars_y = $mkiNy * rnd(_tmp)

 # random variable between 0 and Nstars
 set _tmpvar = $Nstars * rnd(_tmp)

 # random mags drawn from random variable
 set mki_stars_mag = $M0 + log(_tmpvar*$alpha/0.43)/$alpha

 set mki_stars_flux = ten(-0.4*mki_stars_mag)
 sort mki_stars_mag mki_stars_flux mki_stars_x mki_stars_y

 # insert the stars
 mcreate mki_stars_image $mkiNx $mkiNy
 for i 0 $Nstars
   zap mki_stars_image mki_stars_x[$i] mki_stars_y[$i] 1 1 -v mki_stars_flux[$i]
 end
end

####
macro mkimage.smooth
 if ($0 != 5)
   echo "USAGE: smooth (in) (out) (sigma) (func)"
   break
 end
 fft2d $1 0 to sm_I_r sm_I_i
 set sm_g = zero($1)
 keyword sm_g NAXIS1 SMnx
 keyword sm_g NAXIS2 SMny
 $4 sm_g $3 -c 0 0
 $4 sm_g $3 -c $SMnx 0
 $4 sm_g $3 -c 0 $SMny
 $4 sm_g $3 -c $SMnx $SMny
 stats -q sm_g - - - -
 set sm_g = sm_g/$TOTAL
 fft2d sm_g 0 to sm_G_r sm_G_i
 set sm_r = sm_G_r*sm_I_r - sm_G_i*sm_I_i
 set sm_i = sm_G_r*sm_I_i + sm_G_i*sm_I_r
 fft2d -inverse sm_r sm_i to sm_I_r sm_I_i

 set $2 = sm_I_r*$SMnx*$SMny*$SMnx*$SMny

 if (1)
  delete sm_I_r 
  delete sm_I_i
  delete sm_G_r 
  delete sm_G_i
  delete sm_r
  delete sm_i
  delete sm_g
 end
end

macro mkpower
 if ($0 != 6)
  echo "USAGE: mkpower (buffer) (sigma) -c X Y"
  break
 end

 delete -q tmp

 keyword $1 NAXIS1 tmpnx
 keyword $1 NAXIS2 tmpny
 mcreate tmp $tmpnx $tmpny
 set tmpx = xramp(tmp) - $4
 set tmpy = yramp(tmp) - $5

# $sigma_y = $2 * $AxisRatio
# $sigma_xy = dsin(2*$Theta)*(

 set tmpz = tmpx^2 / (2 * $2^2) + tmpy^2 / (2 * $2^2)
 set $1 = $1 + 1 / (1 + tmpz + tmpz^2.5)
end

macro mkwing
 if ($0 != 6)
  echo "USAGE: mkpower (buffer) (sigma) -c X Y"
  break
 end

 delete -q tmp

 keyword $1 NAXIS1 tmpnx
 keyword $1 NAXIS2 tmpny
 mcreate tmp $tmpnx $tmpny
 set tmpx = xramp(tmp) - $4
 set tmpy = yramp(tmp) - $5

 set tmpz = tmpx^2 / (2 * $2^2) + tmpy^2 / (2 * $2^2)
 set $1 = $1 + 1 / (1 + tmpz + tmpz^2)
end
