IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Aug 31, 2013, 6:02:26 AM (13 years ago)
Author:
eugene
Message:

merge changes from EAM branch (update tap script for galaxies; count deblended satstars for log message; scale down EXT window a bit for speed (OK?); enable interactive mode for PCM fitting; allow NONCONVERGE extended fit models to succeed; init EXT model fits with EXT_AND_SKY just in case; extra verbosity and various experiments on the index iterations)

Location:
trunk/psphot
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot

  • trunk/psphot/test/tap_psphot_galaxygrid.pro

    r35769 r36086  
    2828$RefOptions = $RefOptions -Df STARS.DENSITY 10.0
    2929$RefOptions = $RefOptions -Df STARS.SIGMA.LIM 0.5
     30$RefOptions = $RefOptions -nx 3000 -ny 3000
    3031
    3132# options for the simulated images (using the refimage for the stars)
    … …  
    3334$FakeOptions = $FakeOptions -exptime 30.0
    3435$FakeOptions = $FakeOptions -D PSF.MODEL PS_MODEL_GAUSS
     36$FakeOptions = $FakeOptions -nx 3000 -ny 3000
    3537 
    3638# sample alternate options:
    … …  
    6567$FakeConfig = $FakeConfig -Df GALAXY.INDEX.MIN 1.66
    6668$FakeConfig = $FakeConfig -Df GALAXY.INDEX.MAX 1.66
    67 $FakeConfig = $FakeConfig -Di GALAXY.GRID.DX 120
    68 $FakeConfig = $FakeConfig -Di GALAXY.GRID.DY 120
     69$FakeConfig = $FakeConfig -Di GALAXY.GRID.DX 300
     70$FakeConfig = $FakeConfig -Di GALAXY.GRID.DY 300
    6971
    7072list fwhm
    … …  
    7577end
    7678
    77 macro go.grid.set.devexp
     79# generate fake images and run psphot on them
     80macro normtest.mkexp.devexp
    7881  $FakeConfig = -camera SIMTEST
    7982  $FakeConfig = $FakeConfig -recipe PPSIM STACKTEST.RUN
    … …  
    8487  $FakeConfig = $FakeConfig -Db GALAXY.FAKE T                        ; # generate a "realistic" distribution of galaxies
    8588  $FakeConfig = $FakeConfig -Df GALAXY.MAG 17.0
     89  $FakeConfig = $FakeConfig -Df GALAXY.GRID.MAG 14.5
    8690  $FakeConfig = $FakeConfig -Db GALAXY.GRID T                        ; # generate a grid of galaxies (constant mag)
    8791  $FakeConfig = $FakeConfig -Df GALAXY.THETA.MIN 0
    8892  $FakeConfig = $FakeConfig -Df GALAXY.THETA.MAX 180
    89   $FakeConfig = $FakeConfig -Di GALAXY.GRID.DX 120
    90   $FakeConfig = $FakeConfig -Di GALAXY.GRID.DY 120
     93  $FakeConfig = $FakeConfig -Di GALAXY.GRID.DX 300
     94  $FakeConfig = $FakeConfig -Di GALAXY.GRID.DY 300
    9195  $FakeConfig = $FakeConfig -Df GALAXY.INDEX.MIN 1.0
    9296  $FakeConfig = $FakeConfig -Df GALAXY.INDEX.MAX 1.0
     97  $BaseConfig = $FakeConfig
    9398
    9499  # $FakeConfig = $FakeConfig -D GALAXY.MODEL PS_MODEL_GAUSS
    … …  
    103108  # $FakeConfig = $FakeConfig -Df GALAXY.INDEX.MAX 1.66
    104109
     110  mkdir normtest
     111
    105112  $Nseq = 0
    106113  foreach type EXP DEV
    107114    foreach Rmajor 3 10 30
    108115      foreach Aratio 0.25 0.5 1.0
    109         foreach fwhm 0.8 1.0 1.5
     116        foreach fwhm 1.0
     117          $FakeConfig = $BaseConfig
    110118          $FakeConfig = $FakeConfig -Df GALAXY.RMAJOR.MIN $Rmajor
    111119          $FakeConfig = $FakeConfig -Df GALAXY.RMAJOR.MAX $Rmajor
    … …  
    113121          $FakeConfig = $FakeConfig -Df GALAXY.ARATIO.MAX $Aratio
    114122         
    115           sprint name "sample.%02d" $Nseq
     123          sprint name "normtest/test.%02d" $Nseq
    116124          mkexp $name $fwhm $type
    117           fitexp $name $name.fit $type\_CONV
    118125          $Nseq ++
    119126        end
    … …  
    123130end
    124131
    125 macro go.grid.set.sersic
     132# generate fake images and run psphot on them
     133macro grid.mkexp.devexp
     134  if ($0 != 2)
     135    echo "USAGE: grid.mkexp.devexp (dir)"
     136    break
     137  end
     138
     139  mkdir $1
     140
    126141  $FakeConfig = -camera SIMTEST
    127142  $FakeConfig = $FakeConfig -recipe PPSIM STACKTEST.RUN
    … …  
    132147  $FakeConfig = $FakeConfig -Db GALAXY.FAKE T                        ; # generate a "realistic" distribution of galaxies
    133148  $FakeConfig = $FakeConfig -Df GALAXY.MAG 17.0
     149  $FakeConfig = $FakeConfig -Df GALAXY.GRID.MAG 14.5
    134150  $FakeConfig = $FakeConfig -Db GALAXY.GRID T                        ; # generate a grid of galaxies (constant mag)
    135151  $FakeConfig = $FakeConfig -Df GALAXY.THETA.MIN 0
    136152  $FakeConfig = $FakeConfig -Df GALAXY.THETA.MAX 180
    137   $FakeConfig = $FakeConfig -Di GALAXY.GRID.DX 120
    138   $FakeConfig = $FakeConfig -Di GALAXY.GRID.DY 120
     153  $FakeConfig = $FakeConfig -Di GALAXY.GRID.DX 300
     154  $FakeConfig = $FakeConfig -Di GALAXY.GRID.DY 300
     155  $FakeConfig = $FakeConfig -Df GALAXY.INDEX.MIN 1.0
     156  $FakeConfig = $FakeConfig -Df GALAXY.INDEX.MAX 1.0
     157  $BaseConfig = $FakeConfig
     158
     159  $Nseq = 0
     160  foreach type EXP DEV
     161    foreach Rmajor 3 10 30
     162      foreach Aratio 0.25 0.5 1.0
     163        foreach fwhm 0.8 1.0 1.5
     164          $FakeConfig = $BaseConfig
     165          $FakeConfig = $FakeConfig -Df GALAXY.RMAJOR.MIN $Rmajor
     166          $FakeConfig = $FakeConfig -Df GALAXY.RMAJOR.MAX $Rmajor
     167          $FakeConfig = $FakeConfig -Df GALAXY.ARATIO.MIN $Aratio
     168          $FakeConfig = $FakeConfig -Df GALAXY.ARATIO.MAX $Aratio
     169         
     170          sprint name "$1/sample.%02d" $Nseq
     171          mkexp $name $fwhm $type
     172          $Nseq ++
     173        end
     174      end
     175    end
     176  end
     177end
     178
     179# generate fake images and run psphot on them
     180macro grid.fitexp.devexp
     181  if ($0 != 3)
     182    echo "USAGE: grid.fitexp.devexp (srcdir) (outdir)"
     183    break
     184  end
     185
     186   #$Nseq = 0
     187   #foreach type EXP DEV
     188
     189  mkdir $2
     190
     191  $Nseq = 0
     192  foreach type EXP DEV
     193    foreach Rmajor 3 10 30
     194      foreach Aratio 0.25 0.5 1.0
     195        foreach fwhm 0.8 1.0 1.5
     196          sprint srcname "$1/sample.%02d" $Nseq
     197          sprint outname "$2/sample.%02d.fit" $Nseq
     198          fitexp $srcname $outname $type\_CONV
     199          $Nseq ++
     200        end
     201      end
     202    end
     203  end
     204end
     205
     206macro grid.load.devexp
     207  if ($0 != 3)
     208    echo "USAGE: grid.load.devexp (srcdir) (fitdir)"
     209    break
     210  end
     211
     212  delete -q Xin_s Yin_s Min_s Tin_s Rin_s rin_s
     213  delete -q Xot_s Yot_s Mot_s Tot_s Rot_s rot_s
     214  delete -q min_S Min_S
     215
     216  $Nseq = 0
     217  foreach type EXP DEV
     218    foreach Rmajor 3 10 30
     219      foreach Aratio 0.25 0.5 1.0
     220        foreach fwhm 0.8 1.0 1.5
     221          sprint name "sample.%02d" $Nseq
     222          cmf.load.concat $1/$name.dat $2/$name.fit.cmf $type
     223          $Nseq ++
     224        end
     225      end
     226    end
     227  end
     228end
     229
     230macro grid.plots.devexp
     231  if ($0 != 2)
     232    echo "USAGE: grid.plot.devexp (version)"
     233    break
     234  end
     235  # things to examine: theta, Rmajor, AR
     236
     237  # go.grid.check.devexp
     238
     239  # check on theta
     240  set dT = Tot_s - Tin_s
     241  set ARin = rin_s / Rin_s
     242  # lim ARin dT; clear; box; plot ARin dT
     243
     244  subset dTx = dT if (ARin < 0.95)
     245  histogram dTx NdT -8 8 0.1 -range dTi
     246  lim -n 0$1 dTi NdT; clear; box; plot -x 1 dTi NdT
     247  label -x "angle offset (degrees, only non-circular)" -y "number count"
     248  resize 700 320
     249
     250  # check on Rmajor
     251  set dR = Rin_s - Rot_s
     252  # lim Rin_s dR; clear; box; plot Rin_s dR
     253  # lim Rot_s dR; clear; box; plot Rot_s dR
     254  # lim Rin_s dR; clear; box; plot Rin_s dR
     255  create n 0 dR[]
     256  set dRf = dR / Rin_s
     257  lim -n 1$1 n -0.5 0.5; clear; box; plot n dRf
     258  label -x sequence -y "1 - R_out| / R_in|"
     259  resize 700 320
     260
     261  # check on A.Ratio
     262  set ARot = rot_s / Rot_s
     263  # lim ARin ARot; clear; box; plot ARin ARot
     264  set fAR = ARot / ARin
     265  lim -n 2$1 n 0.5 1.5; clear; box; plot n fAR
     266  label -x sequence -y "AR_out| / AR_in|"
     267  resize 700 320
     268
     269  # check on magnitude
     270  set dM = Mot_s - Min_s
     271  lim -n 3$1 n -0.5 0.5; clear; box; plot n dM   
     272  label -x sequence -y "M_out| - M_in|"
     273  resize 700 320
     274end
     275
     276macro grid.mkexp.sersic
     277  if ($0 != 2)
     278    echo "USAGE: grid.mkexp.devexp (dir)"
     279    break
     280  end
     281
     282  mkdir $1
     283
     284  $FakeConfig = -camera SIMTEST
     285  $FakeConfig = $FakeConfig -recipe PPSIM STACKTEST.RUN
     286  $FakeConfig = $FakeConfig -D PSASTRO:PSASTRO.CATDIR catdir.ref
     287  $FakeConfig = $FakeConfig -Db STARS.FAKE F                         ; # only use stars from catdir.ref
     288  $FakeConfig = $FakeConfig -Db MATCH.DENSITY F
     289  $FakeConfig = $FakeConfig -Db PSF.CONVOLVE T
     290  $FakeConfig = $FakeConfig -Db GALAXY.FAKE T                        ; # generate a "realistic" distribution of galaxies
     291  $FakeConfig = $FakeConfig -Df GALAXY.MAG 17.0
     292  $FakeConfig = $FakeConfig -Df GALAXY.GRID.MAG 14.5
     293  $FakeConfig = $FakeConfig -Db GALAXY.GRID T                        ; # generate a grid of galaxies (constant mag)
     294  $FakeConfig = $FakeConfig -Df GALAXY.THETA.MIN 0
     295  $FakeConfig = $FakeConfig -Df GALAXY.THETA.MAX 180
     296  $FakeConfig = $FakeConfig -Di GALAXY.GRID.DX 300
     297  $FakeConfig = $FakeConfig -Di GALAXY.GRID.DY 300
    139298  $FakeConfig = $FakeConfig -D GALAXY.MODEL PS_MODEL_SERSIC
    140 
    141   # $FakeConfig = $FakeConfig -Df GALAXY.INDEX.MIN 1.0
    142   # $FakeConfig = $FakeConfig -Df GALAXY.INDEX.MAX 1.0
    143   # $FakeConfig = $FakeConfig -Df GALAXY.RMAJOR.MIN 10.0
    144   # $FakeConfig = $FakeConfig -Df GALAXY.RMAJOR.MAX 10.0
    145   # $FakeConfig = $FakeConfig -Df GALAXY.ARATIO.MIN 0.25
    146   # $FakeConfig = $FakeConfig -Df GALAXY.ARATIO.MAX 0.25
    147   # $FakeConfig = $FakeConfig -Df GALAXY.INDEX.MIN 1.66
    148   # $FakeConfig = $FakeConfig -Df GALAXY.INDEX.MAX 1.66
     299  $BaseConfig = $FakeConfig
    149300
    150301  $Nseq = 0
    … …  
    153304      foreach Aratio 0.25 0.5 1.0
    154305        foreach fwhm 0.8 1.0 1.5
     306          $FakeConfig = $BaseConfig
    155307          $FakeConfig = $FakeConfig -Df GALAXY.RMAJOR.MIN $Rmajor
    156308          $FakeConfig = $FakeConfig -Df GALAXY.RMAJOR.MAX $Rmajor
    … …  
    160312          $FakeConfig = $FakeConfig -Df GALAXY.INDEX.MAX $index
    161313         
    162           sprint name "sersic.%02d" $Nseq
     314          sprint name "$1/sersic.%02d" $Nseq
    163315          mkexp $name $fwhm SERSIC
    164           fitexp $name $name.fit SER\_CONV
     316          echo "$Nseq : $index $Rmajor $Aratio $fwhm"
    165317          $Nseq ++
    166318        end
    … …  
    170322end
    171323
    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
     324macro grid.fitexp.sersic.devexp
     325  if ($0 != 3)
     326    echo "USAGE: grid.fitexp.sersic.devexp (srcdir) (outdir)"
     327    break
     328  end
     329
     330  mkdir $2
    175331
    176332  $Nseq = 0
    177   foreach type EXP DEV
     333  foreach index 1 2 3 4
    178334    foreach Rmajor 3 10 30
    179335      foreach Aratio 0.25 0.5 1.0
    180336        foreach fwhm 0.8 1.0 1.5
    181           sprint name "sample.%02d" $Nseq
    182           ckgalaxy.load $name.dat $name.fit.cmf
     337          sprint name "sersic.%02d" $Nseq
     338          fitexp $1/$name $2/$name.fit EXP_CONV,DEV_CONV
    183339          $Nseq ++
    184340        end
    … …  
    188344end
    189345
    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
     346macro grid.fitexp.sersic
     347  if ($0 != 3)
     348    echo "USAGE: grid.fitexp.sersic (srcdir) (outdir)"
     349    break
     350  end
     351
     352  mkdir $2
     353
     354  $Nseq = 0
     355  foreach index 1 2 3 4
     356    foreach Rmajor 3 10 30
     357      foreach Aratio 0.25 0.5 1.0
     358        foreach fwhm 0.8 1.0 1.5
     359          sprint name "sersic.%02d" $Nseq
     360          fitexp $1/$name $2/$name.fit SER\_CONV
     361          $Nseq ++
     362        end
     363      end
     364    end
     365  end
     366end
     367
     368macro grid.load.sersic
     369  if ($0 != 3)
     370    echo "USAGE: grid.load.devexp (srcdir) (fitdir)"
     371    break
     372  end
     373
     374  delete -q Xin_s Yin_s Min_s Tin_s Rin_s rin_s MTin_s Iin_s
     375  delete -q Xot_s Yot_s Mot_s Tot_s Rot_s rot_s MTot_s Iot_s
     376
     377  $Nseq = 0
     378  foreach index 1 2 3 4
     379    foreach Rmajor 3 10 30
     380      foreach Aratio 0.25 0.5 1.0
     381        foreach fwhm 0.8 1.0 1.5
     382          sprint name "sersic.%02d" $Nseq
     383          cmf.load.concat $1/$name.dat $2/$name.fit.cmf SERSIC
     384          $Nseq ++
     385        end
     386      end
     387    end
     388  end
     389end
     390
     391# I want to make plots of Iin_s vs Mkron, Mxx, and similar things
     392# this means I need to be able to join Chip.xfit things against Chip.psf
     393# and to join Chip.xfit(DEV) to Chip.xfit(EXP)
     394
     395# I think I need a generic 'JOIN' function
     396
     397macro grid.plots.sersic
     398  if ($0 != 2)
     399    echo "USAGE: grid.plot.devexp (version)"
     400    break
     401  end
     402  # things to examine: theta, Rmajor, AR
     403
     404  # go.grid.check.devexp
     405
     406  # check on theta
     407  set dT = Tot_s - Tin_s
     408  set ARin = rin_s / Rin_s
     409  # lim ARin dT; clear; box; plot ARin dT
     410
     411  subset dTx = dT if (ARin < 0.95)
     412  histogram dTx NdT -8 8 0.1 -range dTi
     413  lim -n 0$1 dTi NdT; clear; box; plot -x 1 dTi NdT
     414  label -x "angle offset (degrees, only non-circular)" -y "number count"
     415  resize 700 320
     416
     417  # check on Rmajor
     418  set dR = Rin_s - Rot_s
     419  # lim Rin_s dR; clear; box; plot Rin_s dR
     420  # lim Rot_s dR; clear; box; plot Rot_s dR
     421  # lim Rin_s dR; clear; box; plot Rin_s dR
     422  create n 0 dR[]
     423  set dRf = dR / Rin_s
     424  lim -n 1$1 n -0.5 0.5; clear; box; plot n dRf
     425  label -x sequence -y "1 - R_out| / R_in|"
     426  resize 700 320
     427
     428  # check on A.Ratio
     429  set ARot = rot_s / Rot_s
     430  # lim ARin ARot; clear; box; plot ARin ARot
     431  set fAR = ARot / ARin
     432  lim -n 2$1 n 0.5 1.5; clear; box; plot n fAR
     433  label -x sequence -y "AR_out| / AR_in|"
     434  resize 700 320
     435
     436  # check on magnitude
     437  set dM = Mot_s - Min_s
     438  lim -n 3$1 n -0.5 0.5; clear; box; plot n dM   
     439  label -x sequence -y "M_out| - M_in|"
     440  resize 700 320
     441
     442  # check on index
     443  set dI = Iot_s - Iin_s
     444  lim -n 4$1 n -0.5 0.5; clear; box; plot n dI   
     445  label -x sequence -y "I_out| - I_in|"
     446  resize 700 320
    218447end
    219448
    … …  
    280509
    281510  # ppImage / psphot on the output
     511  break -auto off
    282512  echo ppImage $ppImageConfig $psphotConfig -file $basename.fits $outname
    283   exec ppImage $ppImageConfig $psphotConfig -file $basename.fits $outname
    284 end
    285 
    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)"
     513  exec echo ppImage $ppImageConfig $psphotConfig -file $basename.fits $outname >& $outname.log
     514  exec ppImage $ppImageConfig $psphotConfig -file $basename.fits $outname >>& $outname.log
     515  break -auto on
     516end
     517
     518macro cmf.load.concat
     519  if ($0 != 4)
     520    echo "USAGE: cmf.load.concat (dat) (cmf) (type)"
    337521    break
    338522  end
    339523
    340524  data $1
    341   read Xin_all 1 Yin_all 2 Type 4 Min_all 5 RmajIn_all 7 RminIn_all 8 ThetaIn_all 9
     525  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
    342526
    343527  subset Xin = Xin_all if (Type == 1)
    344528  subset Yin = Yin_all if (Type == 1)
    345529  subset Min = Min_all if (Type == 1)
     530  subset Fin = Fin_all if (Type == 1)
     531  set min = -2.5*log(Fin)
    346532
    347533  subset Tin_rad = ThetaIn_all if (Type == 1)
    … …  
    351537  subset RminIn = RminIn_all if (Type == 1)
    352538
     539  subset IndexIn = IndexIn_all if (Type == 1)
     540
     541  $TYPE_S = 83
     542  $TYPE_D = 68
     543  $TYPE_E = 69
     544  if ("$3" == "SERSIC")
     545    $InType = $TYPE_S
     546  end
     547  if ("$3" == "DEV")
     548    $InType = $TYPE_D
     549  end
     550  if ("$3" == "EXP")
     551    $InType = $TYPE_E
     552  end
     553
    353554  data $2
    354   read -fits Chip.xfit X_EXT Y_EXT EXT_INST_MAG EXT_WIDTH_MAJ EXT_WIDTH_MIN EXT_THETA
     555
     556  break -auto off
     557  read -fits Chip.xfit X_EXT Y_EXT EXT_INST_MAG EXT_WIDTH_MAJ EXT_WIDTH_MIN EXT_THETA MODEL_TYPE EXT_PAR_07
     558  $reread = not($STATUS)
     559  break -auto on
     560  if ($reread)
     561    read -fits Chip.xfit X_EXT Y_EXT EXT_INST_MAG EXT_WIDTH_MAJ EXT_WIDTH_MIN EXT_THETA MODEL_TYPE
     562    set EXT_PAR_07 = (MODEL_TYPE:9 == 68)*4 + (MODEL_TYPE:9 == 69)
     563  end
     564
    355565  set EXT_THETA_ALT = EXT_THETA * (EXT_THETA >= 0.0) + (EXT_THETA + 3.14159265) * (EXT_THETA < 0.0)
    356566  set EXT_THETA = EXT_THETA_ALT * 180 / 3.14159265
    … …  
    364574  reindex Yin_m = Yin using index2
    365575
     576  set MTin_m = $InType + zero(Xin_m)
     577  reindex MTot_m = MODEL_TYPE:9 using index1
     578
    366579  reindex Mot_m = EXT_INST_MAG using index1
    367580  reindex Tot_m = EXT_THETA using index1
    … …  
    376589  reindex rin_m = RminIn using index2
    377590 
    378   foreach field X Y M T R r
     591  if ("$3" == "SERSIC")
     592    reindex Iot_m = EXT_PAR_07 using index1
     593    reindex Iin_m = IndexIn using index2
     594   $fields = X Y M T R r MT I
     595  else
     596   $fields = X Y M T R r MT
     597  end
     598 
     599  foreach field $fields
    379600    foreach set in ot
    380601      concat $field\$set\_m $field\$set\_s
    381602    end
    382603  end
    383 end
    384 
    385 macro stats.pair
     604
     605  concat min min_S
     606  concat Min Min_S
     607end
     608
     609macro grid.load.sersic.test
    386610  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 
     611    echo "USAGE: grid.load.devexp.test (srcdir) (fitdir)"
     612    break
     613  end
     614
     615  $fields_bt = X Y M T R r I
     616  $fields_ot = Pmag Kmag Amag IDx MTot
     617
     618  foreach field $fields_bt
     619    foreach set in ot
     620      delete -q $field\$set\_dev_s
     621      delete -q $field\$set\_exp_s
     622    end
     623  end
     624  foreach field $fields_ot
     625    delete -q $field\_dev_s
     626    delete -q $field\_exp_s
     627  end
     628
     629  $Nseq = 0
     630  foreach index 1 2 3 4
     631    foreach Rmajor 3 10 30
     632      foreach Aratio 0.25 0.5 1.0
     633        foreach fwhm 0.8 1.0 1.5
     634          sprint name "sersic.%02d" $Nseq
     635          cmf.load.sersic.test $1/$name.dat $2/$name.fit.cmf
     636          $Nseq ++
     637        end
     638      end
     639    end
     640  end
     641end
     642
     643# I have run DEV and EXP against input models of type SERSIC
     644macro cmf.load.sersic.test
     645  if ($0 != 3)
     646    echo "USAGE: cmf.load.sersic.test (dat) (cmf)"
     647    break
     648  end
     649
     650  # input parameters
    515651  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
     652  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
     653
     654  # galaxies only
     655  subset Xin = Xin_all if (Type == 1)
     656  subset Yin = Yin_all if (Type == 1)
     657  subset Min = Min_all if (Type == 1)
     658  subset Fin = Fin_all if (Type == 1)
     659  set min = -2.5*log(Fin)
     660
     661  subset Tin_rad = ThetaIn_all if (Type == 1)
     662  set Tin = Tin_rad * 180 / 3.14159265
     663
     664  subset RmajIn = RmajIn_all if (Type == 1)
     665  subset RminIn = RminIn_all if (Type == 1)
     666
     667  subset IndexIn = IndexIn_all if (Type == 1)
     668
     669  $TYPE_S = 83
     670  $TYPE_D = 68
     671  $TYPE_E = 69
     672
     673  data $2
     674
     675  # load measured values from xfit
     676  read -fits Chip.xfit IPP_IDET X_EXT Y_EXT EXT_INST_MAG EXT_WIDTH_MAJ EXT_WIDTH_MIN EXT_THETA MODEL_TYPE PSF_INST_MAG AP_MAG KRON_MAG
     677  set EXT_PAR_07 = (MODEL_TYPE:9 == $TYPE_D)*4 + (MODEL_TYPE:9 == $TYPE_E)
     678
     679  set EXT_THETA_ALT = EXT_THETA * (EXT_THETA >= 0.0) + (EXT_THETA + 3.14159265) * (EXT_THETA < 0.0)
     680  set EXT_THETA = EXT_THETA_ALT * 180 / 3.14159265
     681  set IPP_IDET_EXT = IPP_IDET
     682 
     683  match2d X_EXT Y_EXT Xin Yin 1.0 -index1 index1 -index2 index2
     684
     685  reindex Xot_m = X_EXT using index1
     686  reindex Yot_m = Y_EXT using index1
     687
     688  reindex Xin_m = Xin using index2
     689  reindex Yin_m = Yin using index2
     690
     691  reindex Mot_m = EXT_INST_MAG using index1
     692  reindex Tot_m = EXT_THETA using index1
     693
     694  reindex Min_m = Min using index2
     695  reindex Tin_m = Tin using index2
     696
     697  reindex Rot_m = EXT_WIDTH_MAJ using index1
     698  reindex rot_m = EXT_WIDTH_MIN using index1
     699
     700  reindex Rin_m = RmajIn using index2
     701  reindex rin_m = RminIn using index2
     702 
     703  reindex Pmag_m = PSF_INST_MAG using index1
     704  reindex Kmag_m = KRON_MAG using index1
     705  reindex Amag_m = AP_MAG using index1
     706
     707  reindex IDx_m = IPP_IDET_EXT using index1
     708
     709  reindex MTot_m = MODEL_TYPE:9 using index1
     710
     711  reindex Iot_m = EXT_PAR_07 using index1
     712  reindex Iin_m = IndexIn using index2
     713 
     714  # load moments and other kron values from Chip.psf
     715  read -fits Chip.psf IPP_IDET X_PSF Y_PSF MOMENTS_XX MOMENTS_XY MOMENTS_YY KRON_FLUX_INNER MOMENTS_R1 MOMENTS_RH
     716  set IPP_IDET_PSF = IPP_IDET
     717
     718  join -outer IPP_IDET_PSF IDx_m
     719  reindex Xp = X_PSF using index2
     720  reindex Yp = Y_PSF using index2
     721
     722  reindex Mxx_m = MOMENTS_XX using index2
     723  reindex Mxy_m = MOMENTS_XY using index2
     724  reindex Myy_m = MOMENTS_YY using index2
     725  reindex Mr1_m = MOMENTS_R1 using index2
     726  reindex Mrh_m = MOMENTS_RH using index2
     727  reindex Kfi_m = KRON_FLUX_INNER using index2
     728  set Kmi_m = -2.5*log(Kfi_m)
     729
     730  $fields_bt = X Y M T R r I
     731  $fields_ot = Pmag Kmag Amag IDx MTot Mxx Mxy Myy Mr1 Mrh Kmi
     732
     733  foreach field $fields_ot
     734    subset $field\_exp_m = $field\_m where (MTot_m == $TYPE_E)
     735    subset $field\_dev_m = $field\_m where (MTot_m == $TYPE_D)
     736  end
     737  foreach field $fields_bt
     738    foreach set in ot
     739      subset $field\$set\_exp_m = $field\$set\_m where (MTot_m == $TYPE_E)
     740      subset $field\$set\_dev_m = $field\$set\_m where (MTot_m == $TYPE_D)
     741    end
     742  end
     743
     744  join IDx_exp_m IDx_dev_m
     745  foreach field $fields_ot
     746    reindex $field\_exp_mr = $field\_exp_m using index1
     747    reindex $field\_dev_mr = $field\_dev_m using index2
     748  end
     749  foreach field $fields_bt
     750    foreach set in ot
     751      reindex $field\$set\_exp_mr = $field\$set\_exp_m using index1
     752      reindex $field\$set\_dev_mr = $field\$set\_dev_m using index2
     753    end
     754  end
     755
     756  # concat
     757  foreach field $fields_bt
     758    foreach set in ot
     759      concat $field\$set\_dev_mr $field\$set\_dev_s
     760      concat $field\$set\_exp_mr $field\$set\_exp_s
     761    end
     762  end
     763  foreach field $fields_ot
     764    concat $field\_dev_mr $field\_dev_s
     765    concat $field\_exp_mr $field\_exp_s
     766  end
     767
     768  concat min min_S
     769  concat Min Min_S
    635770end
    636771
    … …  
    659794end
    660795
    661 macro plot.angles
     796macro sersic.integral
    662797  if ($0 != 2)
    663     echo "USAGE: plot.angles (file.dat)"
     798    echo "sersic.integral (index)"
     799    break
     800  end
     801
     802  $index = $1
     803  create r 0.0 200.0 0.001
     804  set q = r^(1/$index)
     805  set f = exp(-q)
     806  set rf = f*r
     807  integrate r rf r[0] r[-1]
     808end
     809
     810macro sersic.integral.rmax
     811  if ($0 != 4)
     812    echo "sersic.integral.rmax (index) (kappa) (rmax)"
     813    echo "kappa is a guess for kappa"
     814    break
     815  end
     816
     817  local index kappa
     818
     819  $index = $1
     820  $kappa = $2
     821
     822  create r 0.0 $3 0.003
     823  set q = r^(1/$index)
     824  set f = exp(-$kappa*q)
     825  set rf = f*r
     826  integrate r rf r[0] r[-1]
     827end
     828
     829# S(r) = exp(-kappa*(r/Reff)^(1/index))
     830# integrate S(r) r dr [ignores 2pi and change-of-variable factors (Rmaj*Rmin)]
     831macro sersic.integral.reff.rmax
     832  if ($0 != 5)
     833    echo "sersic.integral.rmax (index) (kappa) (reff) (rmax)"
     834    echo "kappa is a guess for kappa"
     835    break
     836  end
     837
     838  local index kappa Reff
     839
     840  $index = $1
     841  $kappa = $2
     842  $Reff = $3
     843
     844  create r 0.0 $4 0.01
     845  set q = (r/$Reff)^(1/$index)
     846  set f = exp(-$kappa*q)
     847  set rf = f*r
     848  integrate r rf r[0] r[-1]
     849end
     850
     851# integrate to r = reff (rho = 1), applying kappa
     852macro sersic.integral.reff
     853  if ($0 != 3)
     854    echo "sersic.integral.reff (index) (kappa)"
     855    break
     856  end
     857
     858  local index kappa
     859
     860  $index = $1
     861  $kappa = $2
     862  create r 0.0 100.0 0.001
     863  set q = r^(1/$index)
     864  set f = exp(-$kappa*q)
     865  set rf = f*r
     866  integrate r rf r[0] r[-1]
     867  $S0 = $sum
     868  integrate r rf r[0] 1.0
     869  $S1 = $sum
     870  echo $S0 $S1 {$S1 / $S0}
     871end
     872
     873macro sersic.integral.find.reff
     874  if ($0 != 3)
     875    echo "sersic.integral.find.reff (index) (kappa)"
     876    echo "kappa is a guess for kappa"
     877    break
     878  end
     879
     880  local kappa index Rmax i
     881  $index = $1
     882  $kappa = $1
     883
     884  $Rmin = 0; $Vmin = 0
     885  $Rmax = 4000
     886  sersic.integral.rmax $index $kappa $Rmax
     887  $Vmax = $sum
     888
     889  $Rtry = 0.5*($Rmin + $Rmax)
     890
     891  while (abs($Rmax - $Rmin) > 0.001)
     892    sersic.integral.rmax $index $kappa $Rtry
     893    $Vtry = $sum
     894
     895    $Rold = $Rtry
     896    if ($Vtry > 0.5*$Vmax)
     897      $Rtry = 0.5*($Rmin + $Rtry)
     898      $Rmax = $Rold
     899    else
     900      $Rtry = 0.5*($Rmax + $Rtry)
     901      $Rmin = $Rold
     902    end
     903    # fprintf "%5.2f %5.2f %5.2f | %5.2f | %6.3f %6.3f | %6.3f" $Rmin $Rold $Rmax $Rtry {$Vtry / $Vmax} $Vmax $Vtry
     904  end
     905  # echo $Rtry
     906  # echo {$kappa * $Rtry ^ (1.0/$index)}
     907  $KappaReal = $kappa * $Rtry ^ (1.0/$index)
     908end
     909
     910macro sersic.integral.index
     911
     912  local i kappa
     913
     914  delete -q sersic_sum sersic_norm
     915  # vlist idx 0.5 0.75 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0
     916  # create idx 0.5 5.1 0.1
     917  vlist idx 1 2 4
     918  for i 0 idx[]
     919    $kappa =  -0.275552 + 1.972625*idx[$i] + 0.003487*idx[$i]*idx[$i]
     920    sersic.integral.rmax idx[$i] $kappa 300
     921    concat $sum sersic_sum
     922
     923    $bn = 1.9992*idx[$i] - 0.3271;
     924    $Io = exp($bn);
     925           
     926    # the integral of a Sersic (supposedly) has an analytical form as follows:
     927    $logGamma = lgamma(2.0*idx[$i]);
     928    $bnFactor = $bn^(2.0*idx[$i]);
     929    $norm = idx[$i] * $Io * exp($logGamma) / $bnFactor;
     930    concat $norm sersic_norm
     931
     932    echo idx[$i] $kappa $sum $logGamma $bn $bnFactor $norm
     933  end
     934  set lsersic_sum = ln(sersic_sum)
     935  break
     936 
     937  $order = 2
     938  delete -q sersic_fit fit_idx
     939
     940  # fit to ranges:
     941  delete lsersic_fit fit_idx
     942  vlist bound 0.0 1.0 2.0 3.0 4.0 5.5
     943  for i 0 {bound[] - 1}
     944    subset idxs = idx if (idx >= bound[$i]) && (idx < bound[$i+1])
     945    subset sums = lsersic_sum if (idx >= bound[$i]) && (idx < bound[$i+1])
     946    fit idxs sums $order
     947    applyfit idxs fits
     948    concat fits lsersic_fit
     949    concat idxs fit_idx
     950  end
     951  set sersic_fit = exp(lsersic_fit)
     952
     953  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
     954  lim -n 2 idx lsersic_sum; clear; box; plot idx lsersic_sum; plot -x 0 -c red fit_idx lsersic_fit
     955end
     956
     957macro sersic.kappa.index
     958
     959  local i
     960
     961  delete -q sersic_kappa
     962  # vlist idx 0.5 0.75 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0
     963  create idx 0.5 5.1 0.1
     964  for i 0 idx[]
     965    sersic.integral.find.reff idx[$i] 10
     966    concat $KappaReal sersic_kappa
     967  end
     968 
     969  # $order = 2
     970  # delete -q sersic_kfit kfit_idx
     971  #
     972  # # fit to ranges:
     973  # vlist bound 0.0 1.0 2.0 3.0 4.0 5.5
     974  # for i 0 {bound[] - 1}
     975  #   subset idxs = idx if (idx >= bound[$i]) && (idx < bound[$i+1])
     976  #   subset sums = lsersic_sum if (idx >= bound[$i]) && (idx < bound[$i+1])
     977  #   fit idxs sums $order
     978  #   applyfit idxs fits
     979  #   concat fits lsersic_fit
     980  #   concat idxs fit_idx
     981  # end
     982  # set sersic_fit = exp(lsersic_fit)
     983
     984  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
     985  # lim -n 2 idx lsersic_sum; clear; box; plot idx lsersic_sum; plot -x 0 -c red fit_idx lsersic_fit
     986
     987  # kappa(index) :
     988  # y = -0.275552 x^0 1.972625 x^1 0.003487 x^2
     989end
     990
     991macro sersic.central.pixel
     992  if ($0 != 3)
     993    echo "USAGE: sersic.central.pixel (index) (Reff)"
     994    break
     995  end
     996
     997  local index Reff
     998
     999  $index = $1
     1000  $Reff = $2
     1001
     1002  $kappa = -0.275552 + 1.972625*$index + 0.003487 * $index^2
     1003  # echo "kappa: $kappa"
     1004
     1005  sersic.norm $index
     1006  sersic.integral.reff.rmax $index $kappa $Reff {10*$Reff}; set sumFull = $sum; echo $sum {$Reff^2*$myNorm} {$sum / ($Reff^2*$myNorm)}
     1007  # sersic.integral.reff.rmax $index $kappa $Reff $Reff     ; # echo $sum {$sum / $sumFull}
     1008  # sersic.integral.reff.rmax $index $kappa $Reff 0.1       ; # echo $sum {$sum / $sumFull}
     1009  # sersic.integral.reff.rmax $index $kappa $Reff 0.2       ; # echo $sum {$sum / $sumFull}
     1010  # sersic.integral.reff.rmax $index $kappa $Reff 0.5       ; # echo $sum {$sum / $sumFull}
     1011  # sersic.integral.reff.rmax $index $kappa $Reff 1.0       ; # echo $sum {$sum / $sumFull}
     1012  # sersic.integral.reff.rmax $index $kappa $Reff 2.0       ; # echo $sum {$sum / $sumFull}
     1013  sersic.integral.reff.rmax $index $kappa $Reff 0.564     ; echo $index $Reff {$sum / $sumFull}
     1014
     1015end
     1016
     1017macro sersic.norm
     1018  if ($0 != 2)
     1019    echo "USAGE: sersic.norm (index)"
     1020    break
     1021  end
     1022
     1023  local index
     1024  $index = $1
     1025
     1026  if (($index >= 0.0) && ($index < 1.0))
     1027      $norm = 0.201545  - 0.950965 * $index - 0.315248 * $index^2
     1028      echo $norm {exp($norm)}
     1029      $myNorm = exp($norm)
     1030      return
     1031  end
     1032
     1033  if (($index >= 1.0) && ($index < 2.0))
     1034      $norm = 0.402084  - 1.357775 * $index - 0.105102 * $index^2
     1035      echo $norm {exp($norm)}
     1036      $myNorm = exp($norm)
     1037      return
     1038  end
     1039
     1040  if (($index >= 2.0) && ($index < 3.0))
     1041      $norm = 0.619093 - 1.591674 * $index - 0.041576 * $index^2
     1042      echo $norm {exp($norm)}
     1043      $myNorm = exp($norm)
     1044      return
     1045  end
     1046
     1047  if (($index >= 3.0) && ($index < 4.0))
     1048      $norm = 0.770263 - 1.696421 * $index - 0.023363 * $index^2
     1049      echo $norm {exp($norm)}
     1050      $myNorm = exp($norm)
     1051      return
     1052  end
     1053
     1054  if (($index >= 4.0) && ($index < 5.5))
     1055      $norm = 0.885891 - 1.755684 * $index - 0.015753 * $index^2
     1056      echo $norm {exp($norm)}
     1057      $myNorm = exp($norm)
     1058      return
     1059  end
     1060end
     1061
     1062macro load.model.im
     1063  if ($0 != 2)
     1064    echo "USAGE: load.model.im (N)"
     1065    break
     1066  end
     1067
     1068  rd obj$1 obj.$1.fits
     1069  rd cnv$1 cnv.$1.fits
     1070  rd var$1 var.$1.fits
     1071  rd msk$1 msk.$1.fits
     1072  for i 1 7
     1073    rd dpar$i.$1 dpar.$i.$1.fits
     1074  end
     1075  set dC$1 = (msk$1 == 0) * (obj$1 - cnv$1)^2 / var$1
     1076  tv dC$1 -0.01 3.0
     1077end
     1078
     1079macro load.normdata
     1080  if ($0 != 2)
     1081    echo "USAGE: load.normdata (file)"
    6641082    break
    6651083  end
    6661084
    6671085  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
     1086  read f 3 t 4 m 5
     1087  set M = -2.5*log(f)
     1088  subset Mg = M if (t == 1)
     1089  subset mg = m if (t == 1)
     1090  set dm = mg - Mg
     1091  set n = ramp(dm)
     1092  lim n dm; clear; box; plot n dm
    6961093end
    6971094
Note: See TracChangeset for help on using the changeset viewer.