Changeset 36375 for trunk/psphot/test/tap_psphot_galaxygrid.pro

Timestamp:

Dec 10, 2013, 2:55:11 PM (13 years ago)

Author:

eugene

Message:

merge changes from eam_branches/ipp-20130904

Location:

trunk

Files:

• . (modified) (1 prop)
• psphot (modified) (1 prop)
• psphot/test/tap_psphot_galaxygrid.pro (modified) (14 diffs)

trunk/psphot/test/tap_psphot_galaxygrid.pro

@@ -30 +30 @@
 $RefOptions = $RefOptions -nx 3000 -ny 3000
 
+if (not($?PSFMODEL))
+  $PSFMODEL = PS1_V1
+end
+
+macro reset.options
 # options for the simulated images (using the refimage for the stars)
 $FakeOptions = $BaseOptions
 $FakeOptions = $FakeOptions -exptime 30.0
-$FakeOptions = $FakeOptions -D PSF.MODEL PS_MODEL_GAUSS
+# $FakeOptions = $FakeOptions -D PSF.MODEL PS_MODEL_GAUSS
+$FakeOptions = $FakeOptions -D PSF.MODEL PS_MODEL_$PSFMODEL
 $FakeOptions = $FakeOptions -nx 3000 -ny 3000
   
@@ -69 +75 @@
 $FakeConfig = $FakeConfig -Di GALAXY.GRID.DX 300
 $FakeConfig = $FakeConfig -Di GALAXY.GRID.DY 300
+end
+
+if (not($?FakeConfig)) reset.options
 
 list fwhm 
@@ -75 +84 @@
  1.2 
  1.5
+end
+
+if (not($?CONVOLVE_NSIGMA)) set CONVOLVE_NSIGMA = 5.0
+
+# generate fake images and run psphot on them
+macro mkexp.devexp.single
+  if ($0 != 5)
+    echo "USAGE: mkexp.devexp.single (basename) (type) (Rmajor) (fwhm)"  
+    break
+  end
+
+  $basename = $1
+  $type = $2
+  $Rmajor = $3
+  $fwhm = $4
+
+  $Aratio = 1.0
+
+  $FakeConfig = -camera SIMTEST
+  $FakeConfig = $FakeConfig -recipe PPSIM STACKTEST.RUN
+  $FakeConfig = $FakeConfig -D PSASTRO:PSASTRO.CATDIR catdir.ref
+  $FakeConfig = $FakeConfig -Db STARS.FAKE F                         ; # only use stars from catdir.ref
+  $FakeConfig = $FakeConfig -Db MATCH.DENSITY F
+  $FakeConfig = $FakeConfig -Db PSF.CONVOLVE T
+  $FakeConfig = $FakeConfig -Db GALAXY.FAKE T                        ; # generate a "realistic" distribution of galaxies
+  $FakeConfig = $FakeConfig -Df GALAXY.MAG 17.0
+  $FakeConfig = $FakeConfig -Df GALAXY.GRID.MAG 14.5
+  $FakeConfig = $FakeConfig -Db GALAXY.GRID T                        ; # generate a grid of galaxies (constant mag)
+  $FakeConfig = $FakeConfig -Df GALAXY.THETA.MIN 0 
+  $FakeConfig = $FakeConfig -Df GALAXY.THETA.MAX 180
+  $FakeConfig = $FakeConfig -Di GALAXY.GRID.DX 300
+  $FakeConfig = $FakeConfig -Di GALAXY.GRID.DY 300
+  $FakeConfig = $FakeConfig -Df GALAXY.INDEX.MIN 1.0
+  $FakeConfig = $FakeConfig -Df GALAXY.INDEX.MAX 1.0
+  $FakeConfig = $FakeConfig -Df CONVOLVE.NSIGMA $CONVOLVE_NSIGMA 
+  $BaseConfig = $FakeConfig
+
+  $FakeConfig = $BaseConfig
+  $FakeConfig = $FakeConfig -Df GALAXY.RMAJOR.MIN $Rmajor
+  $FakeConfig = $FakeConfig -Df GALAXY.RMAJOR.MAX $Rmajor
+  $FakeConfig = $FakeConfig -Df GALAXY.ARATIO.MIN $Aratio
+  $FakeConfig = $FakeConfig -Df GALAXY.ARATIO.MAX $Aratio
+          
+  mkexp $basename $fwhm $type
 end
 
@@ -272 +325 @@
   label -x sequence -y "M_out| - M_in|"
   resize 700 320
+
+  # check on magnitude
+  set dI = Iot_s - Iin_s
+  lim -n 4$1 n -1.0 1.0; clear; box; plot n dI
+  label -x sequence -y "I_out| - I_in|"
+  resize 700 320
+end
+
+macro grid.plot.stars
+  if ($0 != 2)
+    echo "USAGE: grid.plot.stars (version)"
+    break
+  end
+  # things to examine: theta, Rmajor, AR
+
+  # go.grid.check.devexp
+
+  # check on position
+  set dX = Xot_s - int(Xin_s) - 0.5
+  set dY = Yot_s - int(Yin_s) - 0.5
+  set ARin = rin_s / Rin_s
+  set dR = Rin_s - Rot_s
+
+  create n 0 dR[]
+  set dRf = dR / Rin_s
+  lim -n 1$1 n -0.5 0.5; clear; box; plot n dRf
+  label -x sequence -y "1 - R_out| / R_in|"
+  resize 700 320
+
+  # check on A.Ratio
+  set ARot = rot_s / Rot_s
+  # lim ARin ARot; clear; box; plot ARin ARot
+  set fAR = ARot / ARin
+  lim -n 2$1 n 0.5 1.5; clear; box; plot n fAR
+  label -x sequence -y "AR_out| / AR_in|"
+  resize 700 320
+
+  # check on magnitude
+  set dM = Mot_s - Min_s
+  lim -n 3$1 n -0.5 0.5; clear; box; plot n dM    
+  label -x sequence -y "M_out| - M_in|"
+  resize 700 320
+
+  lim -n 4$1 n -1.5 1.5; clear; box; plot n dX
+  label -x sequence -y "X_out| - X_in|"
+  resize 700 320
+
+  lim -n 5$1 n -1.5 1.5; clear; box; plot n dY
+  label -x sequence -y "Y_out| - Y_in|"
+  resize 700 320
 end
 
@@ -487 +590 @@
 end
 
+if (not($?NSIGMA_CONV)) set NSIGMA_CONV = 5.0
+
 # create a realistic distribution of fake stars, GAUSS PSF
 macro fitexp
@@ -507 +612 @@
   $psphotConfig = $psphotConfig -Db PSPHOT:SAVE.RESID T
   $psphotConfig = $psphotConfig -D  PSPHOT:EXTENDED_SOURCE_MODELS_SELECTION $fitModel
+  $psphotConfig = $psphotConfig -D  PSPHOT:PSF_MODEL PS_MODEL_$PSFMODEL
+
+  $psphotConfig = $psphotConfig -Db PSPHOT:PSF.RESIDUALS F
+  $psphotConfig = $psphotConfig -Db PSPHOT:POISSON.ERRORS.PHOT.LMM F
+  $psphotConfig = $psphotConfig -Db PSPHOT:EXTENDED_SOURCE_FITS_POISSON F
+  $psphotConfig = $psphotConfig -Di PSPHOT:EXT_FIT_ITER 15
+  $psphotConfig = $psphotConfig -Df PSPHOT:PSF_FIT_RADIUS_SCALE 3.75
+  $psphotConfig = $psphotConfig -Df PSPHOT:EXT_FIT_NSIGMA_CONV $NSIGMA_CONV
 
   # ppImage / psphot on the output
@@ -516 +629 @@
 end
 
+macro cmf.load.reset
+  $fields = X Y M T R r MT I
+  foreach field $fields
+    foreach set in ot
+      delete -q $field\$set\_s
+    end
+  end
+
+  delete -q min_S Min_S
+end
+
 macro cmf.load.concat
   if ($0 != 4)
-    echo "USAGE: cmf.load.concat (dat) (cmf) (type)"
+    echo "USAGE: cmf.load.concat (dat) (cmf) (inType)"
     break
   end
@@ -525 +649 @@
   read Xin_all 1 Yin_all 2 Fin_all 3 Type 4 Min_all 5 RmajIn_all 7 RminIn_all 8 ThetaIn_all 9 IndexIn_all 10
 
+  $TYPE_S = 83
+  $TYPE_D = 68
+  $TYPE_E = 69
+  if ("$3" == "SERSIC")
+    $InType = $TYPE_S
+  end
+  if ("$3" == "DEV")
+    $InType = $TYPE_D
+    set IndexIn_all = 0.125 + zero(Xin_all)
+  end
+  if ("$3" == "EXP")
+    $InType = $TYPE_E
+    set IndexIn_all = 0.5 + zero(Xin_all)
+  end
+
+  # select only the galaxies 
   subset Xin = Xin_all if (Type == 1)
   subset Yin = Yin_all if (Type == 1)
@@ -539 +679 @@
   subset IndexIn = IndexIn_all if (Type == 1)
 
-  $TYPE_S = 83
-  $TYPE_D = 68
-  $TYPE_E = 69
-  if ("$3" == "SERSIC")
-    $InType = $TYPE_S
-  end
-  if ("$3" == "DEV")
-    $InType = $TYPE_D
-  end
-  if ("$3" == "EXP")
-    $InType = $TYPE_E
-  end
-
   data $2
 
   break -auto off
+  output -err /dev/null
   read -fits Chip.xfit X_EXT Y_EXT EXT_INST_MAG EXT_WIDTH_MAJ EXT_WIDTH_MIN EXT_THETA MODEL_TYPE EXT_PAR_07
   $reread = not($STATUS)
+  output -err stderr
   break -auto on
   if ($reread)
     read -fits Chip.xfit X_EXT Y_EXT EXT_INST_MAG EXT_WIDTH_MAJ EXT_WIDTH_MIN EXT_THETA MODEL_TYPE 
-    set EXT_PAR_07 = (MODEL_TYPE:9 == 68)*4 + (MODEL_TYPE:9 == 69)
+    set EXT_PAR_07 = (MODEL_TYPE:9 == 68)*0.125 + (MODEL_TYPE:9 == 69)*0.5
   end
 
@@ -589 +718 @@
   reindex rin_m = RminIn using index2
   
-  if ("$3" == "SERSIC")
-    reindex Iot_m = EXT_PAR_07 using index1
-    reindex Iin_m = IndexIn using index2
-   $fields = X Y M T R r MT I
-  else
-   $fields = X Y M T R r MT
-  end
+  reindex Iot_m = EXT_PAR_07 using index1
+  reindex Iin_m = IndexIn using index2
+  $fields = X Y M T R r MT I
   
   foreach field $fields
@@ -605 +730 @@
   concat min min_S
   concat Min Min_S
+end
+
+macro cmf.load.stars.concat
+  if ($0 != 3)
+    echo "USAGE: cmf.load.concat (dat) (cmf)"
+    break
+  end
+
+  data $1
+  read Xin_all 1 Yin_all 2 Type 4 Min_all 5 RmajIn_all 7 RminIn_all 8 ThetaIn_all 9 IndexIn_all 10
+
+  subset Xin     = Xin_all     if (Type == 0)
+  subset Yin     = Yin_all     if (Type == 0)
+  subset Min     = Min_all     if (Type == 0)
+  subset RmajIn  = RmajIn_all  if (Type == 0)
+  subset RminIn  = RminIn_all  if (Type == 0)
+  subset IndexIn = IndexIn_all if (Type == 0)
+  subset Tin_rad = ThetaIn_all if (Type == 0)
+  set Tin = Tin_rad * 180 / 3.14159265
+
+  data $2
+
+  break -auto off
+  read -fits Chip.psf X_PSF Y_PSF PSF_INST_MAG PSF_MAJOR PSF_MINOR PSF_THETA
+  set PSF_THETA_ALT = PSF_THETA * (PSF_THETA >= 0.0) + (PSF_THETA + 3.14159265) * (PSF_THETA < 0.0)
+  set PSF_THETA = PSF_THETA_ALT * 180 / 3.14159265
+  
+  match2d X_PSF Y_PSF Xin Yin 1.0 -index1 index1 -index2 index2
+
+  reindex Xot_m = X_PSF using index1
+  reindex Yot_m = Y_PSF using index1
+
+  reindex Xin_m = Xin using index2
+  reindex Yin_m = Yin using index2
+
+  reindex Mot_m = PSF_INST_MAG using index1
+  reindex Tot_m = PSF_THETA using index1
+
+  reindex Min_m = Min using index2
+  reindex Tin_m = Tin using index2
+
+  reindex Rot_m = PSF_MAJOR using index1
+  reindex rot_m = PSF_MINOR using index1
+
+  reindex Rin_m = RmajIn using index2
+  reindex rin_m = RminIn using index2
+  
+  $fields = X Y M T R r
+  foreach field $fields
+    foreach set in ot
+      concat $field\$set\_m $field\$set\_s
+    end
+  end
 end
 
@@ -1078 +1256 @@
 
 macro load.normdata
-  if ($0 != 2)
-    echo "USAGE: load.normdata (file)"
+  if ($0 != 3)
+    echo "USAGE: load.normdata (file) [clear/noclear]"
     break
   end
@@ -1090 +1268 @@
   set dm = mg - Mg
   set n = ramp(dm)
-  lim n dm; clear; box; plot n dm
+  if ("$2" == "clear")
+    lim n dm; clear; box; 
+  end
+  plot n dm
 end
 
@@ -1097 +1278 @@
   exit 0
 end
+
+# note that t1 = original confi
+# t2 = no residuals
+# t3 = const weight PSF fit
+# t4 = const weight galaxy fits
+
+macro run.radius.loop
+  local radius Nrun 
+
+  if (1)
+    $Nrun = 0
+    foreach radius 1.0 1.5 2.0 2.5 3.0 4.0 6.0 8.0
+      mkexp.devexp.single tests.20131120/testrad.ps1v1.exp.$Nrun EXP $radius 1.0
+      fitexp tests.20131120/testrad.ps1v1.exp.$Nrun tests.20131120/testrad.ps1v1.exp.$Nrun.t4b EXP_CONV
+      $Nrun ++
+    end
+  end
+
+  $Nrun = 0
+  foreach radius 1.0 1.5 2.0 2.5 3.0 4.0 6.0 8.0
+    cmf.load.concat tests.20131120/testrad.ps1v1.exp.$Nrun.dat tests.20131120/testrad.ps1v1.exp.$Nrun.t4b.cmf EXP
+    echo -no-return $radius ""
+    check.fit tests.20131120/testrad.ps1v1.exp.$Nrun.dat tests.20131120/testrad.ps1v1.exp.$Nrun.t4b.cmf EXP
+    $Nrun ++
+  end
+
+  grid.plots.devexp a
+end
+
+macro check.radius.loop
+  if ($0 != 2)
+    echo "USAGE: check.radius.loop (ext)"
+    break
+  end
+
+  # cmf.load.reset
+
+  local radius Nrun 
+
+  $Nrun = 0
+  foreach radius 1.0 1.5 2.0 2.5 3.0 4.0 6.0 8.0
+    echo -no-return $radius ""
+    check.fit tests.20131120/testrad.ps1v1.exp.$Nrun.dat tests.20131120/testrad.ps1v1.exp.$Nrun.$1.cmf EXP
+    $Nrun ++
+  end
+
+  grid.plots.devexp a
+end
+
+macro check.fit
+  if ($0 != 4)
+    echo "USAGE: check.fit (dat) (cmf) (type)"
+    break
+  end
+
+  cmf.load.concat $1 $2 $3
+  set dM_m = Mot_m - Min_m
+  vstat -q dM_m
+  $mag_off = $MEDIAN
+  set dR_m = Rot_m - Rin_m
+  vstat -q dR_m
+  $rad_off = $MEDIAN
+  echo $mag_off $rad_off
+end
+
+macro make.circles
+
+  delete Xo Yo Ro
+
+  for ix 0 16
+    for iy $ix 16
+      $r2 = $ix^2 + $iy^2
+      if ($r2 > 15^2) continue
+      concat $ix Xo
+      concat $iy Yo
+      concat $r2 Ro
+    end
+  end
+
+  sort Ro Xo Yo
+  $Rold = -1
+  $N = -1
+  for i 0 Ro[]
+    if (Ro[$i] != $Rold)
+      $N ++
+      $Rold = Ro[$i]
+    end
+    if ((Xo[$i] == 0) && (Yo[$i] == 0))
+      fprintf "// center is 0,0"
+      continue
+    end
+    if (Xo[$i] == 0)
+      fprintf "ADD_AXIS (%3d, %2d)     // r^2 = %3d" $N Yo[$i] Ro[$i]
+      continue
+    end
+    if (Xo[$i] == Yo[$i])
+      fprintf "ADD_DIAG (%3d, %2d)     // r^2 = %3d" $N Xo[$i] Ro[$i]
+      continue
+    end
+
+    fprintf "ADD_RAND (%3d, %2d, %2d) // r^2 = %3d" $N Xo[$i] Yo[$i] Ro[$i]
+ end
+end
+
+# EXP : tests.20131120/test.nsig
+# DEV : tests.20131120/test.dev
+
+# run.convolve.loop tests.20131120/test.nsig tests.20131120/test.nsig   EXP EXP_CONV
+# run.convolve.loop tests.20131120/test.dev  tests.20131120/test.dev    DEV DEV_CONV
+# run.convolve.loop tests.20131120/test.nsig tests.20131120/test.serexp EXP SER_CONV
+# run.convolve.loop tests.20131120/test.dev  tests.20131120/test.serdev DEV SER_CONV
+
+macro run.convolve.loop
+  if ($0 != 5)
+    echo "run.convolve.loop (inName) (fitName) (inType) (fitType)"
+    break
+  end
+
+  local radius
+  $radius = 3.0
+  cmf.load.reset
+
+  # foreach Cin 3 5 7 9 11
+  foreach Cin 11
+    $CONVOLVE_NSIGMA = $Cin
+    sprintf nameIn %s.%02d $1 $Cin
+    mkexp.devexp.single $nameIn $3 $radius 1.0
+    foreach Cot 3 5 7 9 11
+      $NSIGMA_CONV = $Cot
+      sprintf nameOt %s.%02d.%02d $2 $Cin $Cot
+      fitexp $nameIn $nameOt $4
+    end
+  end
+
+  # foreach Cin 3 5 7 9 11
+  foreach Cin 11
+    sprintf nameIn %s.%02d $1 $Cin
+    foreach Cot 3 5 7 9 11
+      sprintf nameOt %s.%02d.%02d $2 $Cin $Cot
+      check.fit $nameIn.dat $nameOt.cmf 
+    end
+  end
+  grid.plots.devexp a
+end
+
+macro check.convolve.loop
+  if ($0 != 4)
+    echo "run.convolve.loop (inName) (fitName) (inType)"
+    break
+  end
+
+  local radius
+  $radius = 3.0
+  cmf.load.reset
+
+  # foreach Cin 3 5 7 9 11
+  foreach Cin 11
+    sprintf nameIn %s.%02d $1 $Cin
+    foreach Cot 3 5 7 9 11
+      sprintf nameOt %s.%02d.%02d $2 $Cin $Cot
+      echo -no-return $Cin $Cot " "
+      check.fit $nameIn.dat $nameOt.cmf $3
+    end
+  end
+  grid.plots.devexp a
+end