IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 35973


Ignore:
Timestamp:
Aug 19, 2013, 6:45:43 AM (13 years ago)
Author:
eugene
Message:

on going work on galaxy modeling

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20130711/psphot/test/tap_psphot_galaxygrid.pro

    r35769 r35973  
    7575end
    7676
    77 macro go.grid.set.devexp
     77# generate fake images and run psphot on them
     78macro normtest.mkexp.devexp
    7879  $FakeConfig = -camera SIMTEST
    7980  $FakeConfig = $FakeConfig -recipe PPSIM STACKTEST.RUN
    … …  
    8485  $FakeConfig = $FakeConfig -Db GALAXY.FAKE T                        ; # generate a "realistic" distribution of galaxies
    8586  $FakeConfig = $FakeConfig -Df GALAXY.MAG 17.0
     87  $FakeConfig = $FakeConfig -Df GALAXY.GRID.MAG 14.5
    8688  $FakeConfig = $FakeConfig -Db GALAXY.GRID T                        ; # generate a grid of galaxies (constant mag)
    8789  $FakeConfig = $FakeConfig -Df GALAXY.THETA.MIN 0
    … …  
    9193  $FakeConfig = $FakeConfig -Df GALAXY.INDEX.MIN 1.0
    9294  $FakeConfig = $FakeConfig -Df GALAXY.INDEX.MAX 1.0
     95  $BaseConfig = $FakeConfig
    9396
    9497  # $FakeConfig = $FakeConfig -D GALAXY.MODEL PS_MODEL_GAUSS
    … …  
    103106  # $FakeConfig = $FakeConfig -Df GALAXY.INDEX.MAX 1.66
    104107
     108  mkdir normtest
     109
     110  $Nseq = 0
     111  foreach type EXP DEV
     112    foreach Rmajor 3 10 30
     113      foreach Aratio 0.25 0.5 1.0
     114        foreach fwhm 1.0
     115          $FakeConfig = $BaseConfig
     116          $FakeConfig = $FakeConfig -Df GALAXY.RMAJOR.MIN $Rmajor
     117          $FakeConfig = $FakeConfig -Df GALAXY.RMAJOR.MAX $Rmajor
     118          $FakeConfig = $FakeConfig -Df GALAXY.ARATIO.MIN $Aratio
     119          $FakeConfig = $FakeConfig -Df GALAXY.ARATIO.MAX $Aratio
     120         
     121          sprint name "normtest/test.%02d" $Nseq
     122          mkexp $name $fwhm $type
     123          $Nseq ++
     124        end
     125      end
     126    end
     127  end
     128end
     129
     130# generate fake images and run psphot on them
     131macro grid.mkexp.devexp
     132  if ($0 != 2)
     133    echo "USAGE: grid.mkexp.devexp (dir)"
     134    break
     135  end
     136
     137  mkdir $1
     138
     139  $FakeConfig = -camera SIMTEST
     140  $FakeConfig = $FakeConfig -recipe PPSIM STACKTEST.RUN
     141  $FakeConfig = $FakeConfig -D PSASTRO:PSASTRO.CATDIR catdir.ref
     142  $FakeConfig = $FakeConfig -Db STARS.FAKE F                         ; # only use stars from catdir.ref
     143  $FakeConfig = $FakeConfig -Db MATCH.DENSITY F
     144  $FakeConfig = $FakeConfig -Db PSF.CONVOLVE T
     145  $FakeConfig = $FakeConfig -Db GALAXY.FAKE T                        ; # generate a "realistic" distribution of galaxies
     146  $FakeConfig = $FakeConfig -Df GALAXY.MAG 17.0
     147  $FakeConfig = $FakeConfig -Df GALAXY.GRID.MAG 14.5
     148  $FakeConfig = $FakeConfig -Db GALAXY.GRID T                        ; # generate a grid of galaxies (constant mag)
     149  $FakeConfig = $FakeConfig -Df GALAXY.THETA.MIN 0
     150  $FakeConfig = $FakeConfig -Df GALAXY.THETA.MAX 180
     151  $FakeConfig = $FakeConfig -Di GALAXY.GRID.DX 120
     152  $FakeConfig = $FakeConfig -Di GALAXY.GRID.DY 120
     153  $FakeConfig = $FakeConfig -Df GALAXY.INDEX.MIN 1.0
     154  $FakeConfig = $FakeConfig -Df GALAXY.INDEX.MAX 1.0
     155  $BaseConfig = $FakeConfig
     156
     157  # $FakeConfig = $FakeConfig -D GALAXY.MODEL PS_MODEL_GAUSS
     158  # $FakeConfig = $FakeConfig -D GALAXY.MODEL PS_MODEL_EXP
     159  # $FakeConfig = $FakeConfig -D GALAXY.MODEL PS_MODEL_SERSIC
     160  # $FakeConfig = $FakeConfig -D GALAXY.MODEL PS_MODEL_DEV
     161  # $FakeConfig = $FakeConfig -Df GALAXY.RMAJOR.MIN 10.0
     162  # $FakeConfig = $FakeConfig -Df GALAXY.RMAJOR.MAX 10.0
     163  # $FakeConfig = $FakeConfig -Df GALAXY.ARATIO.MIN 0.25
     164  # $FakeConfig = $FakeConfig -Df GALAXY.ARATIO.MAX 0.25
     165  # $FakeConfig = $FakeConfig -Df GALAXY.INDEX.MIN 1.66
     166  # $FakeConfig = $FakeConfig -Df GALAXY.INDEX.MAX 1.66
     167
    105168  $Nseq = 0
    106169  foreach type EXP DEV
    … …  
    108171      foreach Aratio 0.25 0.5 1.0
    109172        foreach fwhm 0.8 1.0 1.5
     173          $FakeConfig = $BaseConfig
    110174          $FakeConfig = $FakeConfig -Df GALAXY.RMAJOR.MIN $Rmajor
    111175          $FakeConfig = $FakeConfig -Df GALAXY.RMAJOR.MAX $Rmajor
    … …  
    113177          $FakeConfig = $FakeConfig -Df GALAXY.ARATIO.MAX $Aratio
    114178         
    115           sprint name "sample.%02d" $Nseq
     179          sprint name "$1/sample.%02d" $Nseq
    116180          mkexp $name $fwhm $type
    117           fitexp $name $name.fit $type\_CONV
    118181          $Nseq ++
    119182        end
    … …  
    123186end
    124187
    125 macro go.grid.set.sersic
     188# generate fake images and run psphot on them
     189macro grid.fitexp.devexp
     190  if ($0 != 3)
     191    echo "USAGE: grid.fitexp.devexp (srcdir) (outdir)"
     192    break
     193  end
     194
     195   #$Nseq = 0
     196   #foreach type EXP DEV
     197
     198  mkdir $2
     199
     200  $Nseq = 0
     201  foreach type EXP DEV
     202    foreach Rmajor 3 10 30
     203      foreach Aratio 0.25 0.5 1.0
     204        foreach fwhm 0.8 1.0 1.5
     205          sprint srcname "$1/sample.%02d" $Nseq
     206          sprint outname "$2/sample.%02d.fit" $Nseq
     207          fitexp $srcname $outname $type\_CONV
     208          $Nseq ++
     209        end
     210      end
     211    end
     212  end
     213end
     214
     215macro grid.load.devexp
     216  if ($0 != 3)
     217    echo "USAGE: grid.load.devexp (srcdir) (fitdir)"
     218    break
     219  end
     220
     221  delete -q Xin_s Yin_s Min_s Tin_s Rin_s rin_s
     222  delete -q Xot_s Yot_s Mot_s Tot_s Rot_s rot_s
     223  delete -q min_S Min_S
     224
     225  $Nseq = 0
     226  foreach type EXP DEV
     227    foreach Rmajor 3 10 30
     228      foreach Aratio 0.25 0.5 1.0
     229        foreach fwhm 0.8 1.0 1.5
     230          sprint name "sample.%02d" $Nseq
     231          cmf.load.concat $1/$name.dat $2/$name.fit.cmf DEVEXP
     232          $Nseq ++
     233        end
     234      end
     235    end
     236  end
     237end
     238
     239macro grid.plots.devexp
     240  if ($0 != 2)
     241    echo "USAGE: grid.plot.devexp (version)"
     242    break
     243  end
     244  # things to examine: theta, Rmajor, AR
     245
     246  # go.grid.check.devexp
     247
     248  # check on theta
     249  set dT = Tot_s - Tin_s
     250  set ARin = rin_s / Rin_s
     251  # lim ARin dT; clear; box; plot ARin dT
     252
     253  subset dTx = dT if (ARin < 0.95)
     254  histogram dTx NdT -8 8 0.1 -range dTi
     255  lim -n 0$1 dTi NdT; clear; box; plot -x 1 dTi NdT
     256  label -x "angle offset (degrees, only non-circular)" -y "number count"
     257
     258  # check on Rmajor
     259  set dR = Rin_s - Rot_s
     260  # lim Rin_s dR; clear; box; plot Rin_s dR
     261  # lim Rot_s dR; clear; box; plot Rot_s dR
     262  # lim Rin_s dR; clear; box; plot Rin_s dR
     263  create n 0 dR[]
     264  set dRf = dR / Rin_s
     265  lim -n 1$1 n dRf; clear; box; plot n dRf
     266  label -x sequence -y "1 - R_out| / R_in|"
     267
     268  # check on A.Ratio
     269  set ARot = rot_s / Rot_s
     270  # lim ARin ARot; clear; box; plot ARin ARot
     271  set fAR = ARot / ARin
     272  lim -n 2$1 n fAR; clear; box; plot n fAR
     273  label -x sequence -y "AR_out| / AR_in|"
     274
     275  # check on magnitude
     276  set dM = Mot_s - Min_s
     277  lim -n 3$1 n dM; clear; box; plot n dM   
     278  label -x sequence -y "M_out| - M_in|"
     279end
     280
     281macro grid.mkexp.sersic
    126282  $FakeConfig = -camera SIMTEST
    127283  $FakeConfig = $FakeConfig -recipe PPSIM STACKTEST.RUN
    … …  
    138294  $FakeConfig = $FakeConfig -Di GALAXY.GRID.DY 120
    139295  $FakeConfig = $FakeConfig -D GALAXY.MODEL PS_MODEL_SERSIC
     296  $BaseConfig = $FakeConfig
    140297
    141298  # $FakeConfig = $FakeConfig -Df GALAXY.INDEX.MIN 1.0
    … …  
    153310      foreach Aratio 0.25 0.5 1.0
    154311        foreach fwhm 0.8 1.0 1.5
     312          $FakeConfig = $BaseConfig
    155313          $FakeConfig = $FakeConfig -Df GALAXY.RMAJOR.MIN $Rmajor
    156314          $FakeConfig = $FakeConfig -Df GALAXY.RMAJOR.MAX $Rmajor
    … …  
    162320          sprint name "sersic.%02d" $Nseq
    163321          mkexp $name $fwhm SERSIC
     322          $Nseq ++
     323        end
     324      end
     325    end
     326  end
     327end
     328
     329macro grid.fitexp.sersic
     330  $Nseq = 0
     331  foreach index 1 2 3 4
     332    foreach Rmajor 3 10 30
     333      foreach Aratio 0.25 0.5 1.0
     334        foreach fwhm 0.8 1.0 1.5
     335          sprint name "sersic.%02d" $Nseq
    164336          fitexp $name $name.fit SER\_CONV
    165337          $Nseq ++
    … …  
    170342end
    171343
    172 macro go.grid.check.devexp
    173   delete -q Xin_s Yin_s Min_s Tin_s Rin_s rin_s
    174   delete -q Xot_s Yot_s Mot_s Tot_s Rot_s rot_s
     344macro grid.load.sersic
     345  delete -q Xin_s Yin_s Min_s Tin_s Rin_s rin_s Iin_s
     346  delete -q Xot_s Yot_s Mot_s Tot_s Rot_s rot_s Iot_s
    175347
    176348  $Nseq = 0
    177   foreach type EXP DEV
     349  foreach index 1 2 3 4
    178350    foreach Rmajor 3 10 30
    179351      foreach Aratio 0.25 0.5 1.0
    180352        foreach fwhm 0.8 1.0 1.5
    181           sprint name "sample.%02d" $Nseq
    182           ckgalaxy.load $name.dat $name.fit.cmf
     353          sprint name "sersic.%02d" $Nseq
     354          cmf.load.concat $name.dat $name.fit.cmf SERSIC
    183355          $Nseq ++
    184356        end
    … …  
    188360end
    189361
    190 macro go
    191   mkexp test.exp 1.0 EXP
    192   fitexp test.exp test.exp.fit EXP_CONV
    193 
    194   mkexp test.ser 1.0 SERSIC
    195   fitexp test.ser test.ser.fit SER_CONV
    196 
    197   mkexp test.dev 1.0 DEV
    198   fitexp test.dev test.dev.fit DEV_CONV
    199 
    200   mkexp test.gau 1.0 GAUSS
    201   fitexp test.gau test.gau.fit GAU_CONV
    202 
    203   mkexp test.pg 1.0 PGAUSS
    204   fitexp test.pg test.pg.fit PGA_CONV
    205 
    206   mkexp test.qga 1.0 QGAUSS
    207   fitexp test.qga test.qga.fit QGA_CONV
    208 
    209   mkexp test.p1 1.0 PS1_V1
    210   fitexp test.p1 test.p1.fit PS1_CONV
    211 end
    212 
    213 macro go.ckgalaxy
    214   foreach type exp ser dev gau pg qga p1
    215     ckgalaxy test.$type.dat test.$type.fit.cmf
    216     wait $type
    217   end
     362macro grid.plots.sersic
     363  # things to examine: theta, Rmajor, AR
     364
     365  # go.grid.check.devexp
     366
     367  # check on theta
     368  set dT = Tot_s - Tin_s
     369  set ARin = rin_s / Rin_s
     370  # lim ARin dT; clear; box; plot ARin dT
     371
     372  subset dTx = dT if (ARin < 0.95)
     373  histogram dTx NdT -8 8 0.1 -range dTi
     374  lim -n 0 dTi NdT; clear; box; plot -x 1 dTi NdT
     375  label -x "angle offset (degrees, only non-circular)" -y "number count"
     376
     377  # check on Rmajor
     378  set dR = Rin_s - Rot_s
     379  # lim Rin_s dR; clear; box; plot Rin_s dR
     380  # lim Rot_s dR; clear; box; plot Rot_s dR
     381  # lim Rin_s dR; clear; box; plot Rin_s dR
     382  create n 0 dR[]
     383  set dRf = dR / Rin_s
     384  lim -n 1 n dRf; clear; box; plot n dRf
     385  label -x sequence -y "1 - R_out| / R_in|"
     386
     387  # check on A.Ratio
     388  set ARot = rot_s / Rot_s
     389  # lim ARin ARot; clear; box; plot ARin ARot
     390  set fAR = ARot / ARin
     391  lim -n 2 n fAR; clear; box; plot n fAR
     392  label -x sequence -y "AR_out| / AR_in|"
    218393end
    219394
    … …  
    284459end
    285460
    286 macro ckchip
    287   if ($0 != 5)
    288     echo "USAGE: ckchip (raw) (out) (output) (zpt_off)"
    289     break
    290   end
    291 
    292   load.cmf $1 Chip.psf raw
    293   load.cmf $2 Chip.psf out
    294 
    295   # images generated with convolution will not have the right output positions
    296   set X_raw = int(X_PSF_raw) + 0.5
    297   set Y_raw = int(Y_PSF_raw) + 0.5
    298   set M_raw = PSF_INST_MAG_raw + $4
    299   set K_out = -2.5*log(KRON_FLUX_out)
    300   match2d X_PSF_out Y_PSF_out X_PSF_raw Y_PSF_raw 1.5 -index1 index1 -index2 index2
    301 
    302   local i NX NY nx ny N
    303 
    304   device -n compare
    305   resize 1000 1000
    306 
    307   # plot trends as a function of mag
    308   $NX = 2
    309   $NY = 5
    310   $nx = 0
    311   $ny = 0
    312   $N = 0
    313   clear -s
    314   for i 0 $pairs:n
    315     section a$nx\$ny {$nx/$NX} {$ny/$NY} {1/$NX} {1/$NY}
    316     show.pair $i
    317     $ny ++
    318     if ($ny == $NY)
    319       $ny = 0
    320       $nx ++
    321     end
    322     if ($nx == $NX)
    323       png -name $3.$N.png
    324       clear -s
    325       $nx = 0
    326       $ny = 0
    327       $N ++
    328     end
    329   end
    330 
    331   # plot (input - output) vs mag
    332 end
    333 
    334 macro ckgalaxy.load
    335   if ($0 != 3)
    336     echo "USAGE: ckgalaxy (dat) (cmf)"
     461macro cmf.load.concat
     462  if ($0 != 4)
     463    echo "USAGE: cmf.load.concat (dat) (cmf) (type)"
    337464    break
    338465  end
    339466
    340467  data $1
    341   read Xin_all 1 Yin_all 2 Type 4 Min_all 5 RmajIn_all 7 RminIn_all 8 ThetaIn_all 9
     468  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
    342469
    343470  subset Xin = Xin_all if (Type == 1)
    344471  subset Yin = Yin_all if (Type == 1)
    345472  subset Min = Min_all if (Type == 1)
     473  subset Fin = Fin_all if (Type == 1)
     474  set min = -2.5*log(Fin)
    346475
    347476  subset Tin_rad = ThetaIn_all if (Type == 1)
    … …  
    351480  subset RminIn = RminIn_all if (Type == 1)
    352481
     482  subset IndexIn = IndexIn_all if (Type == 1)
     483
    353484  data $2
    354   read -fits Chip.xfit X_EXT Y_EXT EXT_INST_MAG EXT_WIDTH_MAJ EXT_WIDTH_MIN EXT_THETA
     485  if ("$3" == "SERSIC")
     486    read -fits Chip.xfit X_EXT Y_EXT EXT_INST_MAG EXT_WIDTH_MAJ EXT_WIDTH_MIN EXT_THETA EXT_PAR_07
     487  else
     488    read -fits Chip.xfit X_EXT Y_EXT EXT_INST_MAG EXT_WIDTH_MAJ EXT_WIDTH_MIN EXT_THETA
     489  end
    355490  set EXT_THETA_ALT = EXT_THETA * (EXT_THETA >= 0.0) + (EXT_THETA + 3.14159265) * (EXT_THETA < 0.0)
    356491  set EXT_THETA = EXT_THETA_ALT * 180 / 3.14159265
    … …  
    376511  reindex rin_m = RminIn using index2
    377512 
    378   foreach field X Y M T R r
     513  if ("$3" == "SERSIC")
     514    reindex Iot_m = EXT_PAR_07 using index1
     515    reindex Iin_m = IndexIn using index2
     516   $fields = X Y M T R r I
     517  else
     518   $fields = X Y M T R r
     519  end
     520 
     521  foreach field $fields
    379522    foreach set in ot
    380523      concat $field\$set\_m $field\$set\_s
    381524    end
    382525  end
    383 end
    384 
    385 macro stats.pair
    386   if ($0 != 3)
    387     echo "USAGE: stats.pair (N) (output)"
    388     break
    389   end
    390 
    391   list word -split $spairs:$1
    392   if ($word:n != 8)
    393     echo "invalid pair $1"
    394     break
    395   end
    396 
    397   $Nr = $word:3
    398 
    399   reindex v1 = $word:0 using index1
    400   reindex v2 = $word:1 using index2
    401   reindex rv = $word:2 using index$Nr
    402 
    403   set delta = v1 - v2
    404   subset d1 = delta if ($word:4 < rv) && (rv < $word:5) && (abs(delta) < $word:7)
    405   subset d2 = delta if ($word:5 < rv) && (rv < $word:6) && (abs(delta) < $word:7)
    406 
    407   vstats -q d1 -sigma-clip 3.0
    408   $M1 = $MEAN
    409   $S1 = $SIGMA
    410   vstats -q d2 -sigma-clip 3.0
    411   $M2 = $MEAN
    412   $S2 = $SIGMA
    413 
    414   output $2
    415   fprintf "%-18s  %7.4f %7.4f  %7.4f %7.4f" $word:0  $M1 $S1  $M2 $S2
    416   output stdout
    417 end
    418 
    419 macro show.pair
    420   if ($0 != 2)
    421     echo "USAGE: show.pair (N)"
    422     break
    423   end
    424 
    425   list word -split $pairs:$1
    426   if ($word:n != 7)
    427     echo "invalid pair $1"
    428     break
    429   end
    430 
    431   $Nr = $word:3
    432 
    433   reindex v1 = $word:0 using index1
    434   if ("$word:6" == "V")
    435     reindex v2 = $word:1 using index2
    436   end
    437   if ("$word:6" == "S")
    438     set v2 = $word:1 + zero(index1)
    439   end
    440   reindex rv = $word:2 using index$Nr
    441 
    442   set delta = v1 - v2
    443   if (("$word:4" == "def") || ("$word:5" == "def"))
    444     lim rv delta; box; plot rv delta
    445   else
    446     lim rv $word:4 $word:5; box; plot rv delta
    447   end
    448   label -y '$word:0' -x '$word:2'
    449 end
    450 
    451 # This list is used to compare a pair of vectors (sans error) or a
    452 # vector and an expected (constant) value.  The last field defines a
    453 # vector or constant for the comparison.  It is assumed that the
    454 # vector sets have been loaded and matched with match2d to generate
    455 # index vectors 'index1' and index2'.  The macro 'show.pair' generates
    456 # a plot of the range vector vs (v1 - v2).  The indices for v1, v2 are
    457 # index1 and 2 respectively.  The index for the range vector is defined
    458 # by the integer following that vector.  the y-limits of the plot are
    459 # given by the last two numbers
    460 list pairs
    461   X_PSF_out             X_raw                 PSF_INST_MAG_raw 2 -1.01 1.01 V
    462   Y_PSF_out             Y_raw                 PSF_INST_MAG_raw 2 -1.01 1.01 V
    463   X_PSF_out             X_PSF_raw             PSF_INST_MAG_raw 2 -1.01 1.01 V
    464   Y_PSF_out             Y_PSF_raw             PSF_INST_MAG_raw 2 -1.01 1.01 V
    465   X_PSF_SIG_out         X_PSF_SIG_raw         PSF_INST_MAG_raw 2 -1.01 1.01 V
    466   Y_PSF_SIG_out         Y_PSF_SIG_raw         PSF_INST_MAG_raw 2 -1.01 1.01 V
    467   #PSF_INST_MAG_out      PSF_INST_MAG_raw      PSF_INST_MAG_raw 2 -1.01 1.01 V
    468   PSF_INST_MAG_out      M_raw                 PSF_INST_MAG_raw 2 -1.01 1.01 V
    469   PSF_INST_MAG_SIG_out  PSF_INST_MAG_SIG_raw  PSF_INST_MAG_raw 2 -1.01 1.01 V
    470   #PSF_INST_FLUX_out     PSF_INST_FLUX_raw     PSF_INST_MAG_raw 2   def  def V
    471   #PSF_INST_FLUX_SIG_out PSF_INST_FLUX_SIG_raw PSF_INST_MAG_raw 2   def  def V
    472   AP_MAG_out            M_raw                 PSF_INST_MAG_raw 2 -1.01 1.01 V
    473   AP_MAG_RAW_out        M_raw                 PSF_INST_MAG_raw 2 -1.01 1.01 V
    474   AP_MAG_RADIUS_out     0.0                   PSF_INST_MAG_raw 2 -0.01 20.1 S
    475   SKY_out               0.0                   PSF_INST_MAG_raw 2   def  def S
    476   SKY_SIGMA_out         0.0                   PSF_INST_MAG_raw 2   def  def S
    477   PSF_CHISQ_out         1.0                   PSF_INST_MAG_raw 2   def  def S
    478   CR_NSIGMA_out         0.0                   PSF_INST_MAG_raw 2   def  def S
    479   EXT_NSIGMA_out        0.0                   PSF_INST_MAG_raw 2 -5.01 5.01 S
    480   PSF_MAJOR_out         0.0                   PSF_INST_MAG_raw 2 -0.01 5.01 S
    481   PSF_MINOR_out         0.0                   PSF_INST_MAG_raw 2 -0.01 5.01 S
    482   PSF_THETA_out         0.0                   PSF_INST_MAG_raw 2 -1.61 1.61 S
    483   PSF_QF_out            0.0                   PSF_INST_MAG_raw 2 -0.10 1.10 S
    484   PSF_QF_PERFECT_out    0.0                   PSF_INST_MAG_raw 2 -0.10 1.10 S
    485   PSF_NDOF_out          0.0                   PSF_INST_MAG_raw 2  def  def  S
    486   PSF_NPIX_out          0.0                   PSF_INST_MAG_raw 2  def  def  S
    487   MOMENTS_XX_out        0.0                   PSF_INST_MAG_raw 2 -0.01 3.01 S
    488   MOMENTS_XY_out        0.0                   PSF_INST_MAG_raw 2 -3.01 3.01 S
    489   MOMENTS_YY_out        0.0                   PSF_INST_MAG_raw 2 -0.01 3.01 S
    490   MOMENTS_M3C_out       0.0                   PSF_INST_MAG_raw 2 -3.01 3.01 S
    491   MOMENTS_M3S_out       0.0                   PSF_INST_MAG_raw 2 -3.01 3.01 S
    492   MOMENTS_M4C_out       0.0                   PSF_INST_MAG_raw 2 -2.01 2.01 S
    493   MOMENTS_M4S_out       0.0                   PSF_INST_MAG_raw 2 -2.01 2.01 S
    494   MOMENTS_R1_out        0.0                   PSF_INST_MAG_raw 2 -5.01 5.01 S
    495   MOMENTS_RH_out        0.0                   PSF_INST_MAG_raw 2 -5.01 5.01 S
    496   K_out                 M_raw                 PSF_INST_MAG_raw 2  def  def  V
    497   KRON_FLUX_ERR_out     0.0                   PSF_INST_MAG_raw 2  def  def  S
    498   KRON_FLUX_INNER_out   0.0                   PSF_INST_MAG_raw 2  def  def  S
    499   KRON_FLUX_OUTER_out   0.0                   PSF_INST_MAG_raw 2  def  def  S
    500 # CAL_PSF_MAG          CAL_PSF_MAG          none Mraw 2 -1.01 1.01 V
    501 # CAL_PSF_MAG_SIG      CAL_PSF_MAG_SIG      none Mraw 2 -1.01 1.01 V
    502 # RA_PSF               RA_PSF               none Mraw 2 -1.01 1.01 V
    503 # DEC_PSF              DEC_PSF              none Mraw 2 -1.01 1.01 V
    504 # PEAK_FLUX_AS_MAG     PEAK_FLUX_AS_MAG     none Mraw 2 -1.01 1.01 V
    505 # FLAGS                FLAGS                0.0  Mraw 2 -1.01 1.01 S
    506 # FLAGS2               FLAGS2               0.0  Mraw 2 -1.01 1.01 S
    507 end
    508 
    509 macro load.cmf
    510   if ($0 != 4)
    511    echo "load.cmf (filename) (ext) (label)"
    512    break
    513   end
    514 
    515   data $1
    516 
    517   # create the list of fields to load
    518   delete -q myFields
    519   for i 0 $fields:n
    520     if ($?myFields)
    521       $myFields = $myFields $fields:$i
    522     else
    523       $myFields = $fields:$i
    524     end
    525   end
    526 
    527   read -fits $2 $myFields
    528 
    529   # rename the loaded vectors appending the supplied lable
    530   for i 0 $fields:n
    531     set $fields:$i\_$3 = $fields:$i
    532     delete $fields:$i
    533   end
    534 end
    535 
    536 # this list defines the fields to be loaded from file
    537 list fields
    538   X_PSF             
    539   Y_PSF             
    540   X_PSF_SIG         
    541   Y_PSF_SIG         
    542   PSF_INST_MAG       
    543   PSF_INST_MAG_SIG   
    544   PSF_INST_FLUX     
    545   PSF_INST_FLUX_SIG 
    546   AP_MAG             
    547   AP_MAG_RAW         
    548   AP_MAG_RADIUS     
    549   SKY               
    550   SKY_SIGMA         
    551   PSF_CHISQ         
    552   CR_NSIGMA         
    553   EXT_NSIGMA         
    554   PSF_MAJOR         
    555   PSF_MINOR         
    556   PSF_THETA         
    557   PSF_QF             
    558   PSF_QF_PERFECT     
    559   PSF_NDOF           
    560   PSF_NPIX           
    561   MOMENTS_XX         
    562   MOMENTS_XY         
    563   MOMENTS_YY         
    564   MOMENTS_M3C       
    565   MOMENTS_M3S       
    566   MOMENTS_M4C       
    567   MOMENTS_M4S       
    568   MOMENTS_R1         
    569   MOMENTS_RH         
    570   KRON_FLUX         
    571   KRON_FLUX_ERR     
    572   KRON_FLUX_INNER   
    573   KRON_FLUX_OUTER   
    574 #   CAL_PSF_MAG     
    575 #   CAL_PSF_MAG_SIG 
    576 #   RA_PSF           
    577 #   DEC_PSF         
    578 #   PEAK_FLUX_AS_MAG
    579 #  FLAGS           
    580 #  FLAGS2           
    581 end
    582 
    583 # use these cmf entries to measure average stats of the given pairs
    584 list spairs
    585   X_PSF_out             X_raw                 PSF_INST_MAG_raw 2 -15 -10 -8.5  0.1
    586   Y_PSF_out             Y_raw                 PSF_INST_MAG_raw 2 -15 -10 -8.5  0.1
    587   PSF_INST_MAG_out      M_raw                 PSF_INST_MAG_raw 2 -15 -10 -8.5  0.1
    588   AP_MAG_out            M_raw                 PSF_INST_MAG_raw 2 -15 -10 -8.5  0.1
    589 end
    590 
    591 macro show.dpair
    592   if ($0 != 2)
    593     echo "USAGE: show.dpair (N)"
    594     break
    595   end
    596 
    597   list word -split $dpairs:$1
    598   if ($word:n != 7)
    599     echo "invalid dpair $1"
    600     break
    601   end
    602 
    603   $Nr = $word:4
    604 
    605   reindex v1 = $word:0 using index1
    606   reindex dv = $word:1 using index1
    607   reindex v2 = $word:2 using index2
    608   reindex rv = $word:3 using index$Nr
    609 
    610   set delta = (v1 - v2) / dv
    611   if (("$word:5" == "def") || ("$word:6" == "def"))
    612     lim rv delta; box; plot rv delta
    613   else
    614     lim rv $word:5 $word:6; box; plot rv delta
    615   end
    616   label -y '$word:0' -x '$word:3'
    617 end
    618 
    619 # this list is used to compare a pair of vectors with an error it is
    620 # assumed that the vector sets have been loaded and matched with
    621 # match2d to generate index vectors 'index1' and index2'.  the macro
    622 # show.dpair generates a plot of the range vector vs (v1 - v2) / dv1.
    623 # The indices for v1, dv1, and v2 are index1,1, and 2 respectively.  The
    624 # index for the range vector is defined by the integer following that
    625 # vector.  The y-limits of the plot are given by the last two numbers
    626 # (use 'def') for the full default range of the delta vector
    627 list dpairs
    628   # v1              dv                v2   range
    629   X_PSF             X_PSF_SIG         Xraw Mraw  2 -10.0 10.0
    630   Y_PSF             Y_PSF_SIG         Yraw Mraw  2 -10.0 10.0
    631   PSF_INST_MAG      PSF_INST_MAG_SIG  Mraw Mraw  2 -10.0 10.0
    632   PSF_INST_FLUX     PSF_INST_FLUX_SIG Fraw Mraw  2 -10.0 10.0
    633   AP_MAG            PSF_INST_MAG_SIG  Mraw Mraw  2 -10.0 10.0
    634   AP_MAG_RAW        PSF_INST_MAG_SIG  Mraw Mraw  2 -10.0 10.0
     526
     527  concat min min_S
     528  concat Min Min_S
    635529end
    636530
    … …  
    659553end
    660554
    661 macro plot.angles
     555macro sersic.integral
    662556  if ($0 != 2)
    663     echo "USAGE: plot.angles (file.dat)"
     557    echo "sersic.integral (index)"
     558    break
     559  end
     560
     561  $index = $1
     562  create r 0.0 200.0 0.001
     563  set q = r^(1/$index)
     564  set f = exp(-q)
     565  set rf = f*r
     566  integrate r rf r[0] r[-1]
     567end
     568
     569macro sersic.integral.rmax
     570  if ($0 != 4)
     571    echo "sersic.integral.rmax (index) (kappa) (rmax)"
     572    echo "kappa is a guess for kappa"
     573    break
     574  end
     575
     576  local index kappa
     577
     578  $index = $1
     579  $kappa = $2
     580
     581  create r 0.0 $3 0.003
     582  set q = r^(1/$index)
     583  set f = exp(-$kappa*q)
     584  set rf = f*r
     585  integrate r rf r[0] r[-1]
     586end
     587
     588# S(r) = exp(-kappa*(r/Reff)^(1/index))
     589# integrate S(r) r dr [ignores 2pi and change-of-variable factors (Rmaj*Rmin)]
     590macro sersic.integral.reff.rmax
     591  if ($0 != 5)
     592    echo "sersic.integral.rmax (index) (kappa) (reff) (rmax)"
     593    echo "kappa is a guess for kappa"
     594    break
     595  end
     596
     597  local index kappa Reff
     598
     599  $index = $1
     600  $kappa = $2
     601  $Reff = $3
     602
     603  create r 0.0 $4 0.01
     604  set q = (r/$Reff)^(1/$index)
     605  set f = exp(-$kappa*q)
     606  set rf = f*r
     607  integrate r rf r[0] r[-1]
     608end
     609
     610# integrate to r = reff (rho = 1), applying kappa
     611macro sersic.integral.reff
     612  if ($0 != 3)
     613    echo "sersic.integral.reff (index) (kappa)"
     614    break
     615  end
     616
     617  local index kappa
     618
     619  $index = $1
     620  $kappa = $2
     621  create r 0.0 100.0 0.001
     622  set q = r^(1/$index)
     623  set f = exp(-$kappa*q)
     624  set rf = f*r
     625  integrate r rf r[0] r[-1]
     626  $S0 = $sum
     627  integrate r rf r[0] 1.0
     628  $S1 = $sum
     629  echo $S0 $S1 {$S1 / $S0}
     630end
     631
     632macro sersic.integral.find.reff
     633  if ($0 != 3)
     634    echo "sersic.integral.find.reff (index) (kappa)"
     635    echo "kappa is a guess for kappa"
     636    break
     637  end
     638
     639  local kappa index Rmax i
     640  $index = $1
     641  $kappa = $1
     642
     643  $Rmin = 0; $Vmin = 0
     644  $Rmax = 4000
     645  sersic.integral.rmax $index $kappa $Rmax
     646  $Vmax = $sum
     647
     648  $Rtry = 0.5*($Rmin + $Rmax)
     649
     650  while (abs($Rmax - $Rmin) > 0.001)
     651    sersic.integral.rmax $index $kappa $Rtry
     652    $Vtry = $sum
     653
     654    $Rold = $Rtry
     655    if ($Vtry > 0.5*$Vmax)
     656      $Rtry = 0.5*($Rmin + $Rtry)
     657      $Rmax = $Rold
     658    else
     659      $Rtry = 0.5*($Rmax + $Rtry)
     660      $Rmin = $Rold
     661    end
     662    # fprintf "%5.2f %5.2f %5.2f | %5.2f | %6.3f %6.3f | %6.3f" $Rmin $Rold $Rmax $Rtry {$Vtry / $Vmax} $Vmax $Vtry
     663  end
     664  # echo $Rtry
     665  # echo {$kappa * $Rtry ^ (1.0/$index)}
     666  $KappaReal = $kappa * $Rtry ^ (1.0/$index)
     667end
     668
     669macro sersic.integral.index
     670
     671  local i kappa
     672
     673  delete -q sersic_sum sersic_norm
     674  # vlist idx 0.5 0.75 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0
     675  # create idx 0.5 5.1 0.1
     676  vlist idx 1 2 4
     677  for i 0 idx[]
     678    $kappa =  -0.275552 + 1.972625*idx[$i] + 0.003487*idx[$i]*idx[$i]
     679    sersic.integral.rmax idx[$i] $kappa 300
     680    concat $sum sersic_sum
     681
     682    $bn = 1.9992*idx[$i] - 0.3271;
     683    $Io = exp($bn);
     684           
     685    # the integral of a Sersic (supposedly) has an analytical form as follows:
     686    $logGamma = lgamma(2.0*idx[$i]);
     687    $bnFactor = $bn^(2.0*idx[$i]);
     688    $norm = idx[$i] * $Io * exp($logGamma) / $bnFactor;
     689    concat $norm sersic_norm
     690
     691    echo idx[$i] $kappa $sum $logGamma $bn $bnFactor $norm
     692  end
     693  set lsersic_sum = ln(sersic_sum)
     694  break
     695 
     696  $order = 2
     697  delete -q sersic_fit fit_idx
     698
     699  # fit to ranges:
     700  delete lsersic_fit fit_idx
     701  vlist bound 0.0 1.0 2.0 3.0 4.0 5.5
     702  for i 0 {bound[] - 1}
     703    subset idxs = idx if (idx >= bound[$i]) && (idx < bound[$i+1])
     704    subset sums = lsersic_sum if (idx >= bound[$i]) && (idx < bound[$i+1])
     705    fit idxs sums $order
     706    applyfit idxs fits
     707    concat fits lsersic_fit
     708    concat idxs fit_idx
     709  end
     710  set sersic_fit = exp(lsersic_fit)
     711
     712  lim -n 1 idx sersic_sum; clear; box; plot idx sersic_sum; plot -c blue -pt 7 idx sersic_norm; plot -x 0 -c red fit_idx sersic_fit
     713  lim -n 2 idx lsersic_sum; clear; box; plot idx lsersic_sum; plot -x 0 -c red fit_idx lsersic_fit
     714end
     715
     716macro sersic.kappa.index
     717
     718  local i
     719
     720  delete -q sersic_kappa
     721  # vlist idx 0.5 0.75 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0
     722  create idx 0.5 5.1 0.1
     723  for i 0 idx[]
     724    sersic.integral.find.reff idx[$i] 10
     725    concat $KappaReal sersic_kappa
     726  end
     727 
     728  # $order = 2
     729  # delete -q sersic_kfit kfit_idx
     730  #
     731  # # fit to ranges:
     732  # vlist bound 0.0 1.0 2.0 3.0 4.0 5.5
     733  # for i 0 {bound[] - 1}
     734  #   subset idxs = idx if (idx >= bound[$i]) && (idx < bound[$i+1])
     735  #   subset sums = lsersic_sum if (idx >= bound[$i]) && (idx < bound[$i+1])
     736  #   fit idxs sums $order
     737  #   applyfit idxs fits
     738  #   concat fits lsersic_fit
     739  #   concat idxs fit_idx
     740  # end
     741  # set sersic_fit = exp(lsersic_fit)
     742
     743  lim -n 1 idx sersic_kappa; clear; box; plot idx sersic_kappa; # plot -c blue -pt 7 idx sersic_norm; plot -x 0 -c red fit_idx sersic_fit
     744  # lim -n 2 idx lsersic_sum; clear; box; plot idx lsersic_sum; plot -x 0 -c red fit_idx lsersic_fit
     745
     746  # kappa(index) :
     747  # y = -0.275552 x^0 1.972625 x^1 0.003487 x^2
     748end
     749
     750macro sersic.central.pixel
     751  if ($0 != 3)
     752    echo "USAGE: sersic.central.pixel (index) (Reff)"
     753    break
     754  end
     755
     756  local index Reff
     757
     758  $index = $1
     759  $Reff = $2
     760
     761  $kappa = -0.275552 + 1.972625*$index + 0.003487 * $index^2
     762  # echo "kappa: $kappa"
     763
     764  sersic.norm $index
     765  sersic.integral.reff.rmax $index $kappa $Reff {10*$Reff}; set sumFull = $sum; echo $sum {$Reff^2*$myNorm} {$sum / ($Reff^2*$myNorm)}
     766  # sersic.integral.reff.rmax $index $kappa $Reff $Reff     ; # echo $sum {$sum / $sumFull}
     767  # sersic.integral.reff.rmax $index $kappa $Reff 0.1       ; # echo $sum {$sum / $sumFull}
     768  # sersic.integral.reff.rmax $index $kappa $Reff 0.2       ; # echo $sum {$sum / $sumFull}
     769  # sersic.integral.reff.rmax $index $kappa $Reff 0.5       ; # echo $sum {$sum / $sumFull}
     770  # sersic.integral.reff.rmax $index $kappa $Reff 1.0       ; # echo $sum {$sum / $sumFull}
     771  # sersic.integral.reff.rmax $index $kappa $Reff 2.0       ; # echo $sum {$sum / $sumFull}
     772  sersic.integral.reff.rmax $index $kappa $Reff 0.564     ; echo $index $Reff {$sum / $sumFull}
     773
     774end
     775
     776macro sersic.norm
     777  if ($0 != 2)
     778    echo "USAGE: sersic.norm (index)"
     779    break
     780  end
     781
     782  local index
     783  $index = $1
     784
     785  if (($index >= 0.0) && ($index < 1.0))
     786      $norm = 0.201545  - 0.950965 * $index - 0.315248 * $index^2
     787      echo $norm {exp($norm)}
     788      $myNorm = exp($norm)
     789      return
     790  end
     791
     792  if (($index >= 1.0) && ($index < 2.0))
     793      $norm = 0.402084  - 1.357775 * $index - 0.105102 * $index^2
     794      echo $norm {exp($norm)}
     795      $myNorm = exp($norm)
     796      return
     797  end
     798
     799  if (($index >= 2.0) && ($index < 3.0))
     800      $norm = 0.619093 - 1.591674 * $index - 0.041576 * $index^2
     801      echo $norm {exp($norm)}
     802      $myNorm = exp($norm)
     803      return
     804  end
     805
     806  if (($index >= 3.0) && ($index < 4.0))
     807      $norm = 0.770263 - 1.696421 * $index - 0.023363 * $index^2
     808      echo $norm {exp($norm)}
     809      $myNorm = exp($norm)
     810      return
     811  end
     812
     813  if (($index >= 4.0) && ($index < 5.5))
     814      $norm = 0.885891 - 1.755684 * $index - 0.015753 * $index^2
     815      echo $norm {exp($norm)}
     816      $myNorm = exp($norm)
     817      return
     818  end
     819end
     820
     821macro load.model.im
     822  if ($0 != 2)
     823    echo "USAGE: load.model.im (N)"
     824    break
     825  end
     826
     827  rd obj$1 obj.$1.fits
     828  rd cnv$1 cnv.$1.fits
     829  rd var$1 var.$1.fits
     830  rd msk$1 msk.$1.fits
     831  for i 1 7
     832    rd dpar$i.$1 dpar.$i.$1.fits
     833  end
     834  set dC$1 = (msk$1 == 0) * (obj$1 - cnv$1)^2 / var$1
     835  tv dC$1 -0.01 3.0
     836end
     837
     838macro load.normdata
     839  if ($0 != 2)
     840    echo "USAGE: load.normdata (file)"
    664841    break
    665842  end
    666843
    667844  data $1
    668   read x 1 y 2 type 4 trad 9
    669   set t = trad * 180 / 3.14159265
    670   subset xs = x if (type == 1)
    671   subset ys = y if (type == 1)
    672   subset ts = t if (type == 1)
    673 
    674   set cs = dcos(ts)
    675   set sn = dsin(ts)
    676 
    677   delete xp yp
    678   erase red
    679   for i 0 xs[]
    680     point red LINE xs[$i] ys[$i] {10*cs[$i]} {10*sn[$i]}
    681   end
    682   plot xp yp -pt 100
    683 end
    684 
    685 macro check.flux
    686   if ($0 != 2)
    687     echo "USAGE: plot.angles (file.dat)"
    688     break
    689   end
    690 
    691   data $1
    692   read f 3 type 4
    693   subset F = f if (type == 1)
    694   set M = -2.5*log(F)
    695   vstat M
     845  read f 3 t 4 m 5
     846  set M = -2.5*log(f)
     847  subset Mg = M if (t == 1)
     848  subset mg = m if (t == 1)
     849  set dm = mg - Mg
     850  set n = ramp(dm)
     851  lim n dm; clear; box; plot n dm
    696852end
    697853
Note: See TracChangeset for help on using the changeset viewer.