Changeset 40693 for trunk/doc/release.2015/ps1.analysis/examples/moments.sh

Timestamp:

Apr 13, 2019, 12:48:05 PM (7 years ago)

Author:

eugene

Message:

adding smear experiment

File:

• trunk/doc/release.2015/ps1.analysis/examples/moments.sh (modified) (2 diffs)

trunk/doc/release.2015/ps1.analysis/examples/moments.sh

@@ +1 @@
+
+macro imconvolve
+  if ($0 != 4) 
+    echo "USAGE: imconvolve (image) (kernel) (output)"
+    break
+  end
+
+  # if kernel is not the same size as input, pad
+  if ($1[][0] != $2[][0]) 
+    echo "image sizes do not match"
+    break
+  end
+  if ($1[0][] != $2[0][]) 
+    echo "image sizes do not match"
+    break
+  end
+
+  local Nx Ny
+  $Nx = $1[][0]
+  $Ny = $1[0][]
+
+  fft2d $1 0 to im_R im_I
+  fft2d $2 0 to kr_R kr_I
+
+  $Npix = $Nx*$Ny
+  set ff_R = $Npix*(im_R * kr_R - im_I * kr_I)
+  set ff_I = $Npix*(im_I * kr_R + im_R * kr_I)
+
+  fft2d -inverse ff_R ff_I to ot_R ot_I
+  set tmppow = $Npix*sqrt(ot_R^2 + ot_I^2)
+
+  extract tmppow $3       0       0 {$Nx/2} {$Ny/2} {$Nx/2} {$Ny/2} $Nx $Ny
+  extract tmppow $3 {$Nx/2}       0 {$Nx/2} {$Ny/2}       0 {$Ny/2} $Nx $Ny
+  extract tmppow $3       0 {$Ny/2} {$Nx/2} {$Ny/2} {$Nx/2}       0 $Nx $Ny
+  extract tmppow $3 {$Nx/2} {$Ny/2} {$Nx/2} {$Ny/2}       0       0 $Nx $Ny
+end
 
 macro show.pol
@@ -30 +66 @@
 end
 
+macro show.smear
+  if ($0 != 2) 
+    echo "USAGE: show.smear (sigma)"
+    break
+  end
+
+  mcreate z 512 512
+  set x = xramp(z) - z[][0]/2
+  set y = yramp(z) - z[0][]/2
+
+  set f0 = 10*exp(-0.5*((x/2)^2 + (y/10)^2))
+  set psf_r = exp(-0.5*((x/2)^2 + (y/6)^2))
+  rotate psf_r 30.0
+  extract psf_r psf {psf_r[][0]/2 - 256} {psf_r[0][]/2 - 256} 512 512 0 0 512 512
+  star -q psf {psf[][0]/2} {psf[0][]/2} 128
+  $Mxx_p = ($SXg/2.355)^2
+  $Myy_p = ($SYg/2.355)^2
+  $Mxy_p = $SXYg
+
+  delete -q e1_s e2_s e1_o e2_o theta Mxx_s Mxy_s Myy_s Mxx_o Mxy_o Myy_o
+  for rot 0 360 5
+    set fr = f0
+    rotate fr $rot
+    delete -q frs
+    extract fr frs {fr[][0]/2 - 256} {fr[0][]/2 - 256} 512 512 0 0 512 512
+
+    imconvolve frs psf fs
+
+    star -q fs {fs[][0]/2} {fs[0][]/2} 128
+
+    $Mxx = ($SXg/2.355)^2
+    $Myy = ($SYg/2.355)^2
+    $Mxy = $SXYg
+
+    $T = $Mxx + $Myy
+    $e1 = ($Mxx - $Myy) / $T
+    $e2 = 2*$Mxy / $T
+
+    concat $rot theta
+    concat $e1 e1_s
+    concat $e2 e2_s
+    concat $Mxx Mxx_s
+    concat $Mxy Mxy_s
+    concat $Myy Myy_s
+
+    # fr is the rotate version of f0
+    star -q fr {fr[][0]/2} {fr[0][]/2} 128
+
+    $Mxx = ($SXg/2.355)^2
+    $Myy = ($SYg/2.355)^2
+    $Mxy = $SXYg
+
+    $T = $Mxx + $Myy
+    $e1 = ($Mxx - $Myy) / $T
+    $e2 = 2*$Mxy / $T
+
+    concat $e1 e1_o
+    concat $e2 e2_o
+    concat $Mxx Mxx_o
+    concat $Mxy Mxy_o
+    concat $Myy Myy_o
+  end
+
+  lim -n graph -1.5 1.5 -1.5 1.5; clear; box; plot -pt cir -sz 2 -c blue e1_o e2_o; plot -pt cir -sz 2 -c red e1_s e2_s; 
+end
 
 macro test.rot