IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Apr 21, 2014, 5:42:34 AM (12 years ago)
Author:
eugene
Message:

merge changes from trunk

Location:
branches/eam_branches/ps2-tc3-20130727
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ps2-tc3-20130727

  • branches/eam_branches/ps2-tc3-20130727/psphot

    • Property svn:mergeinfo deleted

  • branches/eam_branches/ps2-tc3-20130727/psphot/test/tap_psphot_galaxygrid.pro

    r35769 r36680  
    2828$RefOptions = $RefOptions -Df STARS.DENSITY 10.0
    2929$RefOptions = $RefOptions -Df STARS.SIGMA.LIM 0.5
    30 
     30$RefOptions = $RefOptions -nx 3000 -ny 3000
     31
     32if (not($?PSFMODEL))
     33  $PSFMODEL = PS1_V1
     34end
     35
     36macro reset.options
    3137# options for the simulated images (using the refimage for the stars)
    3238$FakeOptions = $BaseOptions
    3339$FakeOptions = $FakeOptions -exptime 30.0
    34 $FakeOptions = $FakeOptions -D PSF.MODEL PS_MODEL_GAUSS
     40# $FakeOptions = $FakeOptions -D PSF.MODEL PS_MODEL_GAUSS
     41$FakeOptions = $FakeOptions -D PSF.MODEL PS_MODEL_$PSFMODEL
     42$FakeOptions = $FakeOptions -nx 3000 -ny 3000
    3543 
    3644# sample alternate options:
    … …  
    6573$FakeConfig = $FakeConfig -Df GALAXY.INDEX.MIN 1.66
    6674$FakeConfig = $FakeConfig -Df GALAXY.INDEX.MAX 1.66
    67 $FakeConfig = $FakeConfig -Di GALAXY.GRID.DX 120
    68 $FakeConfig = $FakeConfig -Di GALAXY.GRID.DY 120
     75$FakeConfig = $FakeConfig -Di GALAXY.GRID.DX 300
     76$FakeConfig = $FakeConfig -Di GALAXY.GRID.DY 300
     77end
     78
     79if (not($?FakeConfig)) reset.options
    6980
    7081list fwhm
    … …  
    7586end
    7687
    77 macro go.grid.set.devexp
     88if (not($?CONVOLVE_NSIGMA)) set CONVOLVE_NSIGMA = 5.0
     89
     90# generate fake images and run psphot on them
     91macro mkexp.devexp.single
     92  if ($0 != 5)
     93    echo "USAGE: mkexp.devexp.single (basename) (type) (Rmajor) (fwhm)" 
     94    break
     95  end
     96
     97  $basename = $1
     98  $type = $2
     99  $Rmajor = $3
     100  $fwhm = $4
     101
     102  $Aratio = 1.0
     103
    78104  $FakeConfig = -camera SIMTEST
    79105  $FakeConfig = $FakeConfig -recipe PPSIM STACKTEST.RUN
    … …  
    84110  $FakeConfig = $FakeConfig -Db GALAXY.FAKE T                        ; # generate a "realistic" distribution of galaxies
    85111  $FakeConfig = $FakeConfig -Df GALAXY.MAG 17.0
     112  $FakeConfig = $FakeConfig -Df GALAXY.GRID.MAG 14.5
    86113  $FakeConfig = $FakeConfig -Db GALAXY.GRID T                        ; # generate a grid of galaxies (constant mag)
    87114  $FakeConfig = $FakeConfig -Df GALAXY.THETA.MIN 0
    88115  $FakeConfig = $FakeConfig -Df GALAXY.THETA.MAX 180
    89   $FakeConfig = $FakeConfig -Di GALAXY.GRID.DX 120
    90   $FakeConfig = $FakeConfig -Di GALAXY.GRID.DY 120
     116  $FakeConfig = $FakeConfig -Di GALAXY.GRID.DX 300
     117  $FakeConfig = $FakeConfig -Di GALAXY.GRID.DY 300
    91118  $FakeConfig = $FakeConfig -Df GALAXY.INDEX.MIN 1.0
    92119  $FakeConfig = $FakeConfig -Df GALAXY.INDEX.MAX 1.0
     120  $FakeConfig = $FakeConfig -Df CONVOLVE.NSIGMA $CONVOLVE_NSIGMA
     121  $BaseConfig = $FakeConfig
     122
     123  $FakeConfig = $BaseConfig
     124  $FakeConfig = $FakeConfig -Df GALAXY.RMAJOR.MIN $Rmajor
     125  $FakeConfig = $FakeConfig -Df GALAXY.RMAJOR.MAX $Rmajor
     126  $FakeConfig = $FakeConfig -Df GALAXY.ARATIO.MIN $Aratio
     127  $FakeConfig = $FakeConfig -Df GALAXY.ARATIO.MAX $Aratio
     128         
     129  mkexp $basename $fwhm $type
     130end
     131
     132# generate fake images and run psphot on them
     133macro normtest.mkexp.devexp
     134  $FakeConfig = -camera SIMTEST
     135  $FakeConfig = $FakeConfig -recipe PPSIM STACKTEST.RUN
     136  $FakeConfig = $FakeConfig -D PSASTRO:PSASTRO.CATDIR catdir.ref
     137  $FakeConfig = $FakeConfig -Db STARS.FAKE F                         ; # only use stars from catdir.ref
     138  $FakeConfig = $FakeConfig -Db MATCH.DENSITY F
     139  $FakeConfig = $FakeConfig -Db PSF.CONVOLVE T
     140  $FakeConfig = $FakeConfig -Db GALAXY.FAKE T                        ; # generate a "realistic" distribution of galaxies
     141  $FakeConfig = $FakeConfig -Df GALAXY.MAG 17.0
     142  $FakeConfig = $FakeConfig -Df GALAXY.GRID.MAG 14.5
     143  $FakeConfig = $FakeConfig -Db GALAXY.GRID T                        ; # generate a grid of galaxies (constant mag)
     144  $FakeConfig = $FakeConfig -Df GALAXY.THETA.MIN 0
     145  $FakeConfig = $FakeConfig -Df GALAXY.THETA.MAX 180
     146  $FakeConfig = $FakeConfig -Di GALAXY.GRID.DX 300
     147  $FakeConfig = $FakeConfig -Di GALAXY.GRID.DY 300
     148  $FakeConfig = $FakeConfig -Df GALAXY.INDEX.MIN 1.0
     149  $FakeConfig = $FakeConfig -Df GALAXY.INDEX.MAX 1.0
     150  $BaseConfig = $FakeConfig
    93151
    94152  # $FakeConfig = $FakeConfig -D GALAXY.MODEL PS_MODEL_GAUSS
    … …  
    103161  # $FakeConfig = $FakeConfig -Df GALAXY.INDEX.MAX 1.66
    104162
     163  mkdir normtest
     164
    105165  $Nseq = 0
    106166  foreach type EXP DEV
    107167    foreach Rmajor 3 10 30
    108168      foreach Aratio 0.25 0.5 1.0
    109         foreach fwhm 0.8 1.0 1.5
     169        foreach fwhm 1.0
     170          $FakeConfig = $BaseConfig
    110171          $FakeConfig = $FakeConfig -Df GALAXY.RMAJOR.MIN $Rmajor
    111172          $FakeConfig = $FakeConfig -Df GALAXY.RMAJOR.MAX $Rmajor
    … …  
    113174          $FakeConfig = $FakeConfig -Df GALAXY.ARATIO.MAX $Aratio
    114175         
    115           sprint name "sample.%02d" $Nseq
     176          sprint name "normtest/test.%02d" $Nseq
    116177          mkexp $name $fwhm $type
    117           fitexp $name $name.fit $type\_CONV
    118178          $Nseq ++
    119179        end
    … …  
    123183end
    124184
    125 macro go.grid.set.sersic
     185# generate fake images and run psphot on them
     186macro grid.mkexp.devexp
     187  if ($0 != 2)
     188    echo "USAGE: grid.mkexp.devexp (dir)"
     189    break
     190  end
     191
     192  mkdir $1
     193
    126194  $FakeConfig = -camera SIMTEST
    127195  $FakeConfig = $FakeConfig -recipe PPSIM STACKTEST.RUN
    … …  
    132200  $FakeConfig = $FakeConfig -Db GALAXY.FAKE T                        ; # generate a "realistic" distribution of galaxies
    133201  $FakeConfig = $FakeConfig -Df GALAXY.MAG 17.0
     202  $FakeConfig = $FakeConfig -Df GALAXY.GRID.MAG 14.5
    134203  $FakeConfig = $FakeConfig -Db GALAXY.GRID T                        ; # generate a grid of galaxies (constant mag)
    135204  $FakeConfig = $FakeConfig -Df GALAXY.THETA.MIN 0
    136205  $FakeConfig = $FakeConfig -Df GALAXY.THETA.MAX 180
    137   $FakeConfig = $FakeConfig -Di GALAXY.GRID.DX 120
    138   $FakeConfig = $FakeConfig -Di GALAXY.GRID.DY 120
     206  $FakeConfig = $FakeConfig -Di GALAXY.GRID.DX 300
     207  $FakeConfig = $FakeConfig -Di GALAXY.GRID.DY 300
     208  $FakeConfig = $FakeConfig -Df GALAXY.INDEX.MIN 1.0
     209  $FakeConfig = $FakeConfig -Df GALAXY.INDEX.MAX 1.0
     210  $BaseConfig = $FakeConfig
     211
     212  $Nseq = 0
     213  foreach type EXP DEV
     214    foreach Rmajor 3 10 30
     215      foreach Aratio 0.25 0.5 1.0
     216        foreach fwhm 0.8 1.0 1.5
     217          $FakeConfig = $BaseConfig
     218          $FakeConfig = $FakeConfig -Df GALAXY.RMAJOR.MIN $Rmajor
     219          $FakeConfig = $FakeConfig -Df GALAXY.RMAJOR.MAX $Rmajor
     220          $FakeConfig = $FakeConfig -Df GALAXY.ARATIO.MIN $Aratio
     221          $FakeConfig = $FakeConfig -Df GALAXY.ARATIO.MAX $Aratio
     222         
     223          sprint name "$1/sample.%02d" $Nseq
     224          mkexp $name $fwhm $type
     225          $Nseq ++
     226        end
     227      end
     228    end
     229  end
     230end
     231
     232# generate fake images and run psphot on them
     233macro grid.fitexp.devexp
     234  if ($0 != 3)
     235    echo "USAGE: grid.fitexp.devexp (srcdir) (outdir)"
     236    break
     237  end
     238
     239   #$Nseq = 0
     240   #foreach type EXP DEV
     241
     242  mkdir $2
     243
     244  $Nseq = 0
     245  foreach type EXP DEV
     246    foreach Rmajor 3 10 30
     247      foreach Aratio 0.25 0.5 1.0
     248        foreach fwhm 0.8 1.0 1.5
     249          sprint srcname "$1/sample.%02d" $Nseq
     250          sprint outname "$2/sample.%02d.fit" $Nseq
     251          fitexp $srcname $outname $type\_CONV
     252          $Nseq ++
     253        end
     254      end
     255    end
     256  end
     257end
     258
     259macro grid.load.devexp
     260  if ($0 != 3)
     261    echo "USAGE: grid.load.devexp (srcdir) (fitdir)"
     262    break
     263  end
     264
     265  delete -q Xin_s Yin_s Min_s Tin_s Rin_s rin_s
     266  delete -q Xot_s Yot_s Mot_s Tot_s Rot_s rot_s
     267  delete -q min_S Min_S
     268
     269  $Nseq = 0
     270  foreach type EXP DEV
     271    foreach Rmajor 3 10 30
     272      foreach Aratio 0.25 0.5 1.0
     273        foreach fwhm 0.8 1.0 1.5
     274          sprint name "sample.%02d" $Nseq
     275          cmf.load.concat $1/$name.dat $2/$name.fit.cmf $type
     276          $Nseq ++
     277        end
     278      end
     279    end
     280  end
     281end
     282
     283macro grid.plots.devexp
     284  if ($0 != 2)
     285    echo "USAGE: grid.plot.devexp (version)"
     286    break
     287  end
     288  # things to examine: theta, Rmajor, AR
     289
     290  # go.grid.check.devexp
     291
     292  # check on theta
     293  set dT = Tot_s - Tin_s
     294  set ARin = rin_s / Rin_s
     295  # lim ARin dT; clear; box; plot ARin dT
     296
     297  subset dTx = dT if (ARin < 0.95)
     298  histogram dTx NdT -8 8 0.1 -range dTi
     299  lim -n 0$1 dTi NdT; clear; box; plot -x 1 dTi NdT
     300  label -x "angle offset (degrees, only non-circular)" -y "number count"
     301  resize 700 320
     302
     303  # check on Rmajor
     304  set dR = Rin_s - Rot_s
     305  # lim Rin_s dR; clear; box; plot Rin_s dR
     306  # lim Rot_s dR; clear; box; plot Rot_s dR
     307  # lim Rin_s dR; clear; box; plot Rin_s dR
     308  create n 0 dR[]
     309  set dRf = dR / Rin_s
     310  lim -n 1$1 n -0.5 0.5; clear; box; plot n dRf
     311  label -x sequence -y "1 - R_out| / R_in|"
     312  resize 700 320
     313
     314  # check on A.Ratio
     315  set ARot = rot_s / Rot_s
     316  # lim ARin ARot; clear; box; plot ARin ARot
     317  set fAR = ARot / ARin
     318  lim -n 2$1 n 0.5 1.5; clear; box; plot n fAR
     319  label -x sequence -y "AR_out| / AR_in|"
     320  resize 700 320
     321
     322  # check on magnitude
     323  set dM = Mot_s - Min_s
     324  lim -n 3$1 n -0.5 0.5; clear; box; plot n dM   
     325  label -x sequence -y "M_out| - M_in|"
     326  resize 700 320
     327
     328  # check on magnitude
     329  set dI = Iot_s - Iin_s
     330  lim -n 4$1 n -1.0 1.0; clear; box; plot n dI
     331  label -x sequence -y "I_out| - I_in|"
     332  resize 700 320
     333end
     334
     335macro grid.plot.stars
     336  if ($0 != 2)
     337    echo "USAGE: grid.plot.stars (version)"
     338    break
     339  end
     340  # things to examine: theta, Rmajor, AR
     341
     342  # go.grid.check.devexp
     343
     344  # check on position
     345  set dX = Xot_s - int(Xin_s) - 0.5
     346  set dY = Yot_s - int(Yin_s) - 0.5
     347  set ARin = rin_s / Rin_s
     348  set dR = Rin_s - Rot_s
     349
     350  create n 0 dR[]
     351  set dRf = dR / Rin_s
     352  lim -n 1$1 n -0.5 0.5; clear; box; plot n dRf
     353  label -x sequence -y "1 - R_out| / R_in|"
     354  resize 700 320
     355
     356  # check on A.Ratio
     357  set ARot = rot_s / Rot_s
     358  # lim ARin ARot; clear; box; plot ARin ARot
     359  set fAR = ARot / ARin
     360  lim -n 2$1 n 0.5 1.5; clear; box; plot n fAR
     361  label -x sequence -y "AR_out| / AR_in|"
     362  resize 700 320
     363
     364  # check on magnitude
     365  set dM = Mot_s - Min_s
     366  lim -n 3$1 n -0.5 0.5; clear; box; plot n dM   
     367  label -x sequence -y "M_out| - M_in|"
     368  resize 700 320
     369
     370  lim -n 4$1 n -1.5 1.5; clear; box; plot n dX
     371  label -x sequence -y "X_out| - X_in|"
     372  resize 700 320
     373
     374  lim -n 5$1 n -1.5 1.5; clear; box; plot n dY
     375  label -x sequence -y "Y_out| - Y_in|"
     376  resize 700 320
     377end
     378
     379macro grid.mkexp.sersic
     380  if ($0 != 2)
     381    echo "USAGE: grid.mkexp.devexp (dir)"
     382    break
     383  end
     384
     385  mkdir $1
     386
     387  $FakeConfig = -camera SIMTEST
     388  $FakeConfig = $FakeConfig -recipe PPSIM STACKTEST.RUN
     389  $FakeConfig = $FakeConfig -D PSASTRO:PSASTRO.CATDIR catdir.ref
     390  $FakeConfig = $FakeConfig -Db STARS.FAKE F                         ; # only use stars from catdir.ref
     391  $FakeConfig = $FakeConfig -Db MATCH.DENSITY F
     392  $FakeConfig = $FakeConfig -Db PSF.CONVOLVE T
     393  $FakeConfig = $FakeConfig -Db GALAXY.FAKE T                        ; # generate a "realistic" distribution of galaxies
     394  $FakeConfig = $FakeConfig -Df GALAXY.MAG 17.0
     395  $FakeConfig = $FakeConfig -Df GALAXY.GRID.MAG 14.5
     396  $FakeConfig = $FakeConfig -Db GALAXY.GRID T                        ; # generate a grid of galaxies (constant mag)
     397  $FakeConfig = $FakeConfig -Df GALAXY.THETA.MIN 0
     398  $FakeConfig = $FakeConfig -Df GALAXY.THETA.MAX 180
     399  $FakeConfig = $FakeConfig -Di GALAXY.GRID.DX 300
     400  $FakeConfig = $FakeConfig -Di GALAXY.GRID.DY 300
    139401  $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
     402  $BaseConfig = $FakeConfig
    149403
    150404  $Nseq = 0
    … …  
    153407      foreach Aratio 0.25 0.5 1.0
    154408        foreach fwhm 0.8 1.0 1.5
     409          $FakeConfig = $BaseConfig
    155410          $FakeConfig = $FakeConfig -Df GALAXY.RMAJOR.MIN $Rmajor
    156411          $FakeConfig = $FakeConfig -Df GALAXY.RMAJOR.MAX $Rmajor
    … …  
    160415          $FakeConfig = $FakeConfig -Df GALAXY.INDEX.MAX $index
    161416         
    162           sprint name "sersic.%02d" $Nseq
     417          sprint name "$1/sersic.%02d" $Nseq
    163418          mkexp $name $fwhm SERSIC
    164           fitexp $name $name.fit SER\_CONV
     419          echo "$Nseq : $index $Rmajor $Aratio $fwhm"
    165420          $Nseq ++
    166421        end
    … …  
    170425end
    171426
    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
     427macro grid.fitexp.sersic.devexp
     428  if ($0 != 3)
     429    echo "USAGE: grid.fitexp.sersic.devexp (srcdir) (outdir)"
     430    break
     431  end
     432
     433  mkdir $2
    175434
    176435  $Nseq = 0
    177   foreach type EXP DEV
     436  foreach index 1 2 3 4
    178437    foreach Rmajor 3 10 30
    179438      foreach Aratio 0.25 0.5 1.0
    180439        foreach fwhm 0.8 1.0 1.5
    181           sprint name "sample.%02d" $Nseq
    182           ckgalaxy.load $name.dat $name.fit.cmf
     440          sprint name "sersic.%02d" $Nseq
     441          fitexp $1/$name $2/$name.fit EXP_CONV,DEV_CONV
    183442          $Nseq ++
    184443        end
    … …  
    188447end
    189448
    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
     449macro grid.fitexp.sersic
     450  if ($0 != 3)
     451    echo "USAGE: grid.fitexp.sersic (srcdir) (outdir)"
     452    break
     453  end
     454
     455  mkdir $2
     456
     457  $Nseq = 0
     458  foreach index 1 2 3 4
     459    foreach Rmajor 3 10 30
     460      foreach Aratio 0.25 0.5 1.0
     461        foreach fwhm 0.8 1.0 1.5
     462          sprint name "sersic.%02d" $Nseq
     463          fitexp $1/$name $2/$name.fit SER\_CONV
     464          $Nseq ++
     465        end
     466      end
     467    end
     468  end
     469end
     470
     471macro grid.load.sersic
     472  if ($0 != 3)
     473    echo "USAGE: grid.load.devexp (srcdir) (fitdir)"
     474    break
     475  end
     476
     477  delete -q Xin_s Yin_s Min_s Tin_s Rin_s rin_s MTin_s Iin_s
     478  delete -q Xot_s Yot_s Mot_s Tot_s Rot_s rot_s MTot_s Iot_s
     479
     480  $Nseq = 0
     481  foreach index 1 2 3 4
     482    foreach Rmajor 3 10 30
     483      foreach Aratio 0.25 0.5 1.0
     484        foreach fwhm 0.8 1.0 1.5
     485          sprint name "sersic.%02d" $Nseq
     486          cmf.load.concat $1/$name.dat $2/$name.fit.cmf SERSIC
     487          $Nseq ++
     488        end
     489      end
     490    end
     491  end
     492end
     493
     494# I want to make plots of Iin_s vs Mkron, Mxx, and similar things
     495# this means I need to be able to join Chip.xfit things against Chip.psf
     496# and to join Chip.xfit(DEV) to Chip.xfit(EXP)
     497
     498# I think I need a generic 'JOIN' function
     499
     500macro grid.plots.sersic
     501  if ($0 != 2)
     502    echo "USAGE: grid.plot.devexp (version)"
     503    break
     504  end
     505  # things to examine: theta, Rmajor, AR
     506
     507  # go.grid.check.devexp
     508
     509  # check on theta
     510  set dT = Tot_s - Tin_s
     511  set ARin = rin_s / Rin_s
     512  # lim ARin dT; clear; box; plot ARin dT
     513
     514  subset dTx = dT if (ARin < 0.95)
     515  histogram dTx NdT -8 8 0.1 -range dTi
     516  lim -n 0$1 dTi NdT; clear; box; plot -x 1 dTi NdT
     517  label -x "angle offset (degrees, only non-circular)" -y "number count"
     518  resize 700 320
     519
     520  # check on Rmajor
     521  set dR = Rin_s - Rot_s
     522  # lim Rin_s dR; clear; box; plot Rin_s dR
     523  # lim Rot_s dR; clear; box; plot Rot_s dR
     524  # lim Rin_s dR; clear; box; plot Rin_s dR
     525  create n 0 dR[]
     526  set dRf = dR / Rin_s
     527  lim -n 1$1 n -0.5 0.5; clear; box; plot n dRf
     528  label -x sequence -y "1 - R_out| / R_in|"
     529  resize 700 320
     530
     531  # check on A.Ratio
     532  set ARot = rot_s / Rot_s
     533  # lim ARin ARot; clear; box; plot ARin ARot
     534  set fAR = ARot / ARin
     535  lim -n 2$1 n 0.5 1.5; clear; box; plot n fAR
     536  label -x sequence -y "AR_out| / AR_in|"
     537  resize 700 320
     538
     539  # check on magnitude
     540  set dM = Mot_s - Min_s
     541  lim -n 3$1 n -0.5 0.5; clear; box; plot n dM   
     542  label -x sequence -y "M_out| - M_in|"
     543  resize 700 320
     544
     545  # check on index
     546  set dI = Iot_s - Iin_s
     547  lim -n 4$1 n -0.5 0.5; clear; box; plot n dI   
     548  label -x sequence -y "I_out| - I_in|"
     549  resize 700 320
    218550end
    219551
    … …  
    258590end
    259591
     592if (not($?NSIGMA_CONV)) set NSIGMA_CONV = 5.0
     593
    260594# create a realistic distribution of fake stars, GAUSS PSF
    261595macro fitexp
    … …  
    278612  $psphotConfig = $psphotConfig -Db PSPHOT:SAVE.RESID T
    279613  $psphotConfig = $psphotConfig -D  PSPHOT:EXTENDED_SOURCE_MODELS_SELECTION $fitModel
     614  $psphotConfig = $psphotConfig -D  PSPHOT:PSF_MODEL PS_MODEL_$PSFMODEL
     615
     616  $psphotConfig = $psphotConfig -Db PSPHOT:PSF.RESIDUALS F
     617  $psphotConfig = $psphotConfig -Db PSPHOT:POISSON.ERRORS.PHOT.LMM F
     618  $psphotConfig = $psphotConfig -Db PSPHOT:EXTENDED_SOURCE_FITS_POISSON F
     619  $psphotConfig = $psphotConfig -Di PSPHOT:EXT_FIT_ITER 15
     620  $psphotConfig = $psphotConfig -Df PSPHOT:PSF_FIT_RADIUS_SCALE 3.75
     621  $psphotConfig = $psphotConfig -Df PSPHOT:EXT_FIT_NSIGMA_CONV $NSIGMA_CONV
    280622
    281623  # ppImage / psphot on the output
     624  break -auto off
    282625  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)"
     626  exec echo ppImage $ppImageConfig $psphotConfig -file $basename.fits $outname >& $outname.log
     627  exec ppImage $ppImageConfig $psphotConfig -file $basename.fits $outname >>& $outname.log
     628  break -auto on
     629end
     630
     631macro cmf.load.reset
     632  $fields = X Y M T R r MT I
     633  foreach field $fields
     634    foreach set in ot
     635      delete -q $field\$set\_s
     636    end
     637  end
     638
     639  delete -q min_S Min_S
     640end
     641
     642macro cmf.load.concat
     643  if ($0 != 4)
     644    echo "USAGE: cmf.load.concat (dat) (cmf) (inType)"
    337645    break
    338646  end
    339647
    340648  data $1
    341   read Xin_all 1 Yin_all 2 Type 4 Min_all 5 RmajIn_all 7 RminIn_all 8 ThetaIn_all 9
    342 
     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
     650
     651  $TYPE_S = 83
     652  $TYPE_D = 68
     653  $TYPE_E = 69
     654  if ("$3" == "SERSIC")
     655    $InType = $TYPE_S
     656  end
     657  if ("$3" == "DEV")
     658    $InType = $TYPE_D
     659    set IndexIn_all = 0.125 + zero(Xin_all)
     660  end
     661  if ("$3" == "EXP")
     662    $InType = $TYPE_E
     663    set IndexIn_all = 0.5 + zero(Xin_all)
     664  end
     665
     666  # select only the galaxies
    343667  subset Xin = Xin_all if (Type == 1)
    344668  subset Yin = Yin_all if (Type == 1)
    345669  subset Min = Min_all if (Type == 1)
     670  subset Fin = Fin_all if (Type == 1)
     671  set min = -2.5*log(Fin)
    346672
    347673  subset Tin_rad = ThetaIn_all if (Type == 1)
    … …  
    351677  subset RminIn = RminIn_all if (Type == 1)
    352678
     679  subset IndexIn = IndexIn_all if (Type == 1)
     680
    353681  data $2
    354   read -fits Chip.xfit X_EXT Y_EXT EXT_INST_MAG EXT_WIDTH_MAJ EXT_WIDTH_MIN EXT_THETA
     682
     683  break -auto off
     684  output -err /dev/null
     685  read -fits Chip.xfit X_EXT Y_EXT EXT_INST_MAG EXT_WIDTH_MAJ EXT_WIDTH_MIN EXT_THETA MODEL_TYPE EXT_PAR_07
     686  $reread = not($STATUS)
     687  output -err stderr
     688  break -auto on
     689  if ($reread)
     690    read -fits Chip.xfit X_EXT Y_EXT EXT_INST_MAG EXT_WIDTH_MAJ EXT_WIDTH_MIN EXT_THETA MODEL_TYPE
     691    set EXT_PAR_07 = (MODEL_TYPE:9 == 68)*0.125 + (MODEL_TYPE:9 == 69)*0.5
     692  end
     693
    355694  set EXT_THETA_ALT = EXT_THETA * (EXT_THETA >= 0.0) + (EXT_THETA + 3.14159265) * (EXT_THETA < 0.0)
    356695  set EXT_THETA = EXT_THETA_ALT * 180 / 3.14159265
    … …  
    364703  reindex Yin_m = Yin using index2
    365704
     705  set MTin_m = $InType + zero(Xin_m)
     706  reindex MTot_m = MODEL_TYPE:9 using index1
     707
    366708  reindex Mot_m = EXT_INST_MAG using index1
    367709  reindex Tot_m = EXT_THETA using index1
    … …  
    376718  reindex rin_m = RminIn using index2
    377719 
    378   foreach field X Y M T R r
     720  reindex Iot_m = EXT_PAR_07 using index1
     721  reindex Iin_m = IndexIn using index2
     722  $fields = X Y M T R r MT I
     723 
     724  foreach field $fields
    379725    foreach set in ot
    380726      concat $field\$set\_m $field\$set\_s
    381727    end
    382728  end
    383 end
    384 
    385 macro stats.pair
     729
     730  concat min min_S
     731  concat Min Min_S
     732end
     733
     734macro cmf.load.stars.concat
    386735  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
     736    echo "USAGE: cmf.load.concat (dat) (cmf)"
     737    break
    513738  end
    514739
    515740  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
     741  read Xin_all 1 Yin_all 2 Type 4 Min_all 5 RmajIn_all 7 RminIn_all 8 ThetaIn_all 9 IndexIn_all 10
     742
     743  subset Xin     = Xin_all     if (Type == 0)
     744  subset Yin     = Yin_all     if (Type == 0)
     745  subset Min     = Min_all     if (Type == 0)
     746  subset RmajIn  = RmajIn_all  if (Type == 0)
     747  subset RminIn  = RminIn_all  if (Type == 0)
     748  subset IndexIn = IndexIn_all if (Type == 0)
     749  subset Tin_rad = ThetaIn_all if (Type == 0)
     750  set Tin = Tin_rad * 180 / 3.14159265
     751
     752  data $2
     753
     754  break -auto off
     755  read -fits Chip.psf X_PSF Y_PSF PSF_INST_MAG PSF_MAJOR PSF_MINOR PSF_THETA
     756  set PSF_THETA_ALT = PSF_THETA * (PSF_THETA >= 0.0) + (PSF_THETA + 3.14159265) * (PSF_THETA < 0.0)
     757  set PSF_THETA = PSF_THETA_ALT * 180 / 3.14159265
     758 
     759  match2d X_PSF Y_PSF Xin Yin 1.0 -index1 index1 -index2 index2
     760
     761  reindex Xot_m = X_PSF using index1
     762  reindex Yot_m = Y_PSF using index1
     763
     764  reindex Xin_m = Xin using index2
     765  reindex Yin_m = Yin using index2
     766
     767  reindex Mot_m = PSF_INST_MAG using index1
     768  reindex Tot_m = PSF_THETA using index1
     769
     770  reindex Min_m = Min using index2
     771  reindex Tin_m = Tin using index2
     772
     773  reindex Rot_m = PSF_MAJOR using index1
     774  reindex rot_m = PSF_MINOR using index1
     775
     776  reindex Rin_m = RmajIn using index2
     777  reindex rin_m = RminIn using index2
     778 
     779  $fields = X Y M T R r
     780  foreach field $fields
     781    foreach set in ot
     782      concat $field\$set\_m $field\$set\_s
     783    end
     784  end
     785end
     786
     787macro grid.load.sersic.test
     788  if ($0 != 3)
     789    echo "USAGE: grid.load.devexp.test (srcdir) (fitdir)"
     790    break
     791  end
     792
     793  $fields_bt = X Y M T R r I
     794  $fields_ot = Pmag Kmag Amag IDx MTot
     795
     796  foreach field $fields_bt
     797    foreach set in ot
     798      delete -q $field\$set\_dev_s
     799      delete -q $field\$set\_exp_s
     800    end
     801  end
     802  foreach field $fields_ot
     803    delete -q $field\_dev_s
     804    delete -q $field\_exp_s
     805  end
     806
     807  $Nseq = 0
     808  foreach index 1 2 3 4
     809    foreach Rmajor 3 10 30
     810      foreach Aratio 0.25 0.5 1.0
     811        foreach fwhm 0.8 1.0 1.5
     812          sprint name "sersic.%02d" $Nseq
     813          cmf.load.sersic.test $1/$name.dat $2/$name.fit.cmf
     814          $Nseq ++
     815        end
     816      end
     817    end
     818  end
     819end
     820
     821# I have run DEV and EXP against input models of type SERSIC
     822macro cmf.load.sersic.test
     823  if ($0 != 3)
     824    echo "USAGE: cmf.load.sersic.test (dat) (cmf)"
     825    break
     826  end
     827
     828  # input parameters
     829  data $1
     830  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
     831
     832  # galaxies only
     833  subset Xin = Xin_all if (Type == 1)
     834  subset Yin = Yin_all if (Type == 1)
     835  subset Min = Min_all if (Type == 1)
     836  subset Fin = Fin_all if (Type == 1)
     837  set min = -2.5*log(Fin)
     838
     839  subset Tin_rad = ThetaIn_all if (Type == 1)
     840  set Tin = Tin_rad * 180 / 3.14159265
     841
     842  subset RmajIn = RmajIn_all if (Type == 1)
     843  subset RminIn = RminIn_all if (Type == 1)
     844
     845  subset IndexIn = IndexIn_all if (Type == 1)
     846
     847  $TYPE_S = 83
     848  $TYPE_D = 68
     849  $TYPE_E = 69
     850
     851  data $2
     852
     853  # load measured values from xfit
     854  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
     855  set EXT_PAR_07 = (MODEL_TYPE:9 == $TYPE_D)*4 + (MODEL_TYPE:9 == $TYPE_E)
     856
     857  set EXT_THETA_ALT = EXT_THETA * (EXT_THETA >= 0.0) + (EXT_THETA + 3.14159265) * (EXT_THETA < 0.0)
     858  set EXT_THETA = EXT_THETA_ALT * 180 / 3.14159265
     859  set IPP_IDET_EXT = IPP_IDET
     860 
     861  match2d X_EXT Y_EXT Xin Yin 1.0 -index1 index1 -index2 index2
     862
     863  reindex Xot_m = X_EXT using index1
     864  reindex Yot_m = Y_EXT using index1
     865
     866  reindex Xin_m = Xin using index2
     867  reindex Yin_m = Yin using index2
     868
     869  reindex Mot_m = EXT_INST_MAG using index1
     870  reindex Tot_m = EXT_THETA using index1
     871
     872  reindex Min_m = Min using index2
     873  reindex Tin_m = Tin using index2
     874
     875  reindex Rot_m = EXT_WIDTH_MAJ using index1
     876  reindex rot_m = EXT_WIDTH_MIN using index1
     877
     878  reindex Rin_m = RmajIn using index2
     879  reindex rin_m = RminIn using index2
     880 
     881  reindex Pmag_m = PSF_INST_MAG using index1
     882  reindex Kmag_m = KRON_MAG using index1
     883  reindex Amag_m = AP_MAG using index1
     884
     885  reindex IDx_m = IPP_IDET_EXT using index1
     886
     887  reindex MTot_m = MODEL_TYPE:9 using index1
     888
     889  reindex Iot_m = EXT_PAR_07 using index1
     890  reindex Iin_m = IndexIn using index2
     891 
     892  # load moments and other kron values from Chip.psf
     893  read -fits Chip.psf IPP_IDET X_PSF Y_PSF MOMENTS_XX MOMENTS_XY MOMENTS_YY KRON_FLUX_INNER MOMENTS_R1 MOMENTS_RH
     894  set IPP_IDET_PSF = IPP_IDET
     895
     896  join -outer IPP_IDET_PSF IDx_m
     897  reindex Xp = X_PSF using index2
     898  reindex Yp = Y_PSF using index2
     899
     900  reindex Mxx_m = MOMENTS_XX using index2
     901  reindex Mxy_m = MOMENTS_XY using index2
     902  reindex Myy_m = MOMENTS_YY using index2
     903  reindex Mr1_m = MOMENTS_R1 using index2
     904  reindex Mrh_m = MOMENTS_RH using index2
     905  reindex Kfi_m = KRON_FLUX_INNER using index2
     906  set Kmi_m = -2.5*log(Kfi_m)
     907
     908  $fields_bt = X Y M T R r I
     909  $fields_ot = Pmag Kmag Amag IDx MTot Mxx Mxy Myy Mr1 Mrh Kmi
     910
     911  foreach field $fields_ot
     912    subset $field\_exp_m = $field\_m where (MTot_m == $TYPE_E)
     913    subset $field\_dev_m = $field\_m where (MTot_m == $TYPE_D)
     914  end
     915  foreach field $fields_bt
     916    foreach set in ot
     917      subset $field\$set\_exp_m = $field\$set\_m where (MTot_m == $TYPE_E)
     918      subset $field\$set\_dev_m = $field\$set\_m where (MTot_m == $TYPE_D)
     919    end
     920  end
     921
     922  join IDx_exp_m IDx_dev_m
     923  foreach field $fields_ot
     924    reindex $field\_exp_mr = $field\_exp_m using index1
     925    reindex $field\_dev_mr = $field\_dev_m using index2
     926  end
     927  foreach field $fields_bt
     928    foreach set in ot
     929      reindex $field\$set\_exp_mr = $field\$set\_exp_m using index1
     930      reindex $field\$set\_dev_mr = $field\$set\_dev_m using index2
     931    end
     932  end
     933
     934  # concat
     935  foreach field $fields_bt
     936    foreach set in ot
     937      concat $field\$set\_dev_mr $field\$set\_dev_s
     938      concat $field\$set\_exp_mr $field\$set\_exp_s
     939    end
     940  end
     941  foreach field $fields_ot
     942    concat $field\_dev_mr $field\_dev_s
     943    concat $field\_exp_mr $field\_exp_s
     944  end
     945
     946  concat min min_S
     947  concat Min Min_S
    635948end
    636949
    … …  
    659972end
    660973
    661 macro plot.angles
     974macro sersic.integral
    662975  if ($0 != 2)
    663     echo "USAGE: plot.angles (file.dat)"
     976    echo "sersic.integral (index)"
     977    break
     978  end
     979
     980  $index = $1
     981  create r 0.0 200.0 0.001
     982  set q = r^(1/$index)
     983  set f = exp(-q)
     984  set rf = f*r
     985  integrate r rf r[0] r[-1]
     986end
     987
     988macro sersic.integral.rmax
     989  if ($0 != 4)
     990    echo "sersic.integral.rmax (index) (kappa) (rmax)"
     991    echo "kappa is a guess for kappa"
     992    break
     993  end
     994
     995  local index kappa
     996
     997  $index = $1
     998  $kappa = $2
     999
     1000  create r 0.0 $3 0.003
     1001  set q = r^(1/$index)
     1002  set f = exp(-$kappa*q)
     1003  set rf = f*r
     1004  integrate r rf r[0] r[-1]
     1005end
     1006
     1007# S(r) = exp(-kappa*(r/Reff)^(1/index))
     1008# integrate S(r) r dr [ignores 2pi and change-of-variable factors (Rmaj*Rmin)]
     1009macro sersic.integral.reff.rmax
     1010  if ($0 != 5)
     1011    echo "sersic.integral.rmax (index) (kappa) (reff) (rmax)"
     1012    echo "kappa is a guess for kappa"
     1013    break
     1014  end
     1015
     1016  local index kappa Reff
     1017
     1018  $index = $1
     1019  $kappa = $2
     1020  $Reff = $3
     1021
     1022  create r 0.0 $4 0.01
     1023  set q = (r/$Reff)^(1/$index)
     1024  set f = exp(-$kappa*q)
     1025  set rf = f*r
     1026  integrate r rf r[0] r[-1]
     1027end
     1028
     1029# integrate to r = reff (rho = 1), applying kappa
     1030macro sersic.integral.reff
     1031  if ($0 != 3)
     1032    echo "sersic.integral.reff (index) (kappa)"
     1033    break
     1034  end
     1035
     1036  local index kappa
     1037
     1038  $index = $1
     1039  $kappa = $2
     1040  create r 0.0 100.0 0.001
     1041  set q = r^(1/$index)
     1042  set f = exp(-$kappa*q)
     1043  set rf = f*r
     1044  integrate r rf r[0] r[-1]
     1045  $S0 = $sum
     1046  integrate r rf r[0] 1.0
     1047  $S1 = $sum
     1048  echo $S0 $S1 {$S1 / $S0}
     1049end
     1050
     1051macro sersic.integral.find.reff
     1052  if ($0 != 3)
     1053    echo "sersic.integral.find.reff (index) (kappa)"
     1054    echo "kappa is a guess for kappa"
     1055    break
     1056  end
     1057
     1058  local kappa index Rmax i
     1059  $index = $1
     1060  $kappa = $1
     1061
     1062  $Rmin = 0; $Vmin = 0
     1063  $Rmax = 4000
     1064  sersic.integral.rmax $index $kappa $Rmax
     1065  $Vmax = $sum
     1066
     1067  $Rtry = 0.5*($Rmin + $Rmax)
     1068
     1069  while (abs($Rmax - $Rmin) > 0.001)
     1070    sersic.integral.rmax $index $kappa $Rtry
     1071    $Vtry = $sum
     1072
     1073    $Rold = $Rtry
     1074    if ($Vtry > 0.5*$Vmax)
     1075      $Rtry = 0.5*($Rmin + $Rtry)
     1076      $Rmax = $Rold
     1077    else
     1078      $Rtry = 0.5*($Rmax + $Rtry)
     1079      $Rmin = $Rold
     1080    end
     1081    # fprintf "%5.2f %5.2f %5.2f | %5.2f | %6.3f %6.3f | %6.3f" $Rmin $Rold $Rmax $Rtry {$Vtry / $Vmax} $Vmax $Vtry
     1082  end
     1083  # echo $Rtry
     1084  # echo {$kappa * $Rtry ^ (1.0/$index)}
     1085  $KappaReal = $kappa * $Rtry ^ (1.0/$index)
     1086end
     1087
     1088macro sersic.integral.index
     1089
     1090  local i kappa
     1091
     1092  delete -q sersic_sum sersic_norm
     1093  # vlist idx 0.5 0.75 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0
     1094  # create idx 0.5 5.1 0.1
     1095  vlist idx 1 2 4
     1096  for i 0 idx[]
     1097    $kappa =  -0.275552 + 1.972625*idx[$i] + 0.003487*idx[$i]*idx[$i]
     1098    sersic.integral.rmax idx[$i] $kappa 300
     1099    concat $sum sersic_sum
     1100
     1101    $bn = 1.9992*idx[$i] - 0.3271;
     1102    $Io = exp($bn);
     1103           
     1104    # the integral of a Sersic (supposedly) has an analytical form as follows:
     1105    $logGamma = lgamma(2.0*idx[$i]);
     1106    $bnFactor = $bn^(2.0*idx[$i]);
     1107    $norm = idx[$i] * $Io * exp($logGamma) / $bnFactor;
     1108    concat $norm sersic_norm
     1109
     1110    echo idx[$i] $kappa $sum $logGamma $bn $bnFactor $norm
     1111  end
     1112  set lsersic_sum = ln(sersic_sum)
     1113  break
     1114 
     1115  $order = 2
     1116  delete -q sersic_fit fit_idx
     1117
     1118  # fit to ranges:
     1119  delete lsersic_fit fit_idx
     1120  vlist bound 0.0 1.0 2.0 3.0 4.0 5.5
     1121  for i 0 {bound[] - 1}
     1122    subset idxs = idx if (idx >= bound[$i]) && (idx < bound[$i+1])
     1123    subset sums = lsersic_sum if (idx >= bound[$i]) && (idx < bound[$i+1])
     1124    fit idxs sums $order
     1125    applyfit idxs fits
     1126    concat fits lsersic_fit
     1127    concat idxs fit_idx
     1128  end
     1129  set sersic_fit = exp(lsersic_fit)
     1130
     1131  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
     1132  lim -n 2 idx lsersic_sum; clear; box; plot idx lsersic_sum; plot -x 0 -c red fit_idx lsersic_fit
     1133end
     1134
     1135macro sersic.kappa.index
     1136
     1137  local i
     1138
     1139  delete -q sersic_kappa
     1140  # vlist idx 0.5 0.75 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0
     1141  create idx 0.5 5.1 0.1
     1142  for i 0 idx[]
     1143    sersic.integral.find.reff idx[$i] 10
     1144    concat $KappaReal sersic_kappa
     1145  end
     1146 
     1147  # $order = 2
     1148  # delete -q sersic_kfit kfit_idx
     1149  #
     1150  # # fit to ranges:
     1151  # vlist bound 0.0 1.0 2.0 3.0 4.0 5.5
     1152  # for i 0 {bound[] - 1}
     1153  #   subset idxs = idx if (idx >= bound[$i]) && (idx < bound[$i+1])
     1154  #   subset sums = lsersic_sum if (idx >= bound[$i]) && (idx < bound[$i+1])
     1155  #   fit idxs sums $order
     1156  #   applyfit idxs fits
     1157  #   concat fits lsersic_fit
     1158  #   concat idxs fit_idx
     1159  # end
     1160  # set sersic_fit = exp(lsersic_fit)
     1161
     1162  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
     1163  # lim -n 2 idx lsersic_sum; clear; box; plot idx lsersic_sum; plot -x 0 -c red fit_idx lsersic_fit
     1164
     1165  # kappa(index) :
     1166  # y = -0.275552 x^0 1.972625 x^1 0.003487 x^2
     1167end
     1168
     1169macro sersic.central.pixel
     1170  if ($0 != 3)
     1171    echo "USAGE: sersic.central.pixel (index) (Reff)"
     1172    break
     1173  end
     1174
     1175  local index Reff
     1176
     1177  $index = $1
     1178  $Reff = $2
     1179
     1180  $kappa = -0.275552 + 1.972625*$index + 0.003487 * $index^2
     1181  # echo "kappa: $kappa"
     1182
     1183  sersic.norm $index
     1184  sersic.integral.reff.rmax $index $kappa $Reff {10*$Reff}; set sumFull = $sum; echo $sum {$Reff^2*$myNorm} {$sum / ($Reff^2*$myNorm)}
     1185  # sersic.integral.reff.rmax $index $kappa $Reff $Reff     ; # echo $sum {$sum / $sumFull}
     1186  # sersic.integral.reff.rmax $index $kappa $Reff 0.1       ; # echo $sum {$sum / $sumFull}
     1187  # sersic.integral.reff.rmax $index $kappa $Reff 0.2       ; # echo $sum {$sum / $sumFull}
     1188  # sersic.integral.reff.rmax $index $kappa $Reff 0.5       ; # echo $sum {$sum / $sumFull}
     1189  # sersic.integral.reff.rmax $index $kappa $Reff 1.0       ; # echo $sum {$sum / $sumFull}
     1190  # sersic.integral.reff.rmax $index $kappa $Reff 2.0       ; # echo $sum {$sum / $sumFull}
     1191  sersic.integral.reff.rmax $index $kappa $Reff 0.564     ; echo $index $Reff {$sum / $sumFull}
     1192
     1193end
     1194
     1195macro sersic.norm
     1196  if ($0 != 2)
     1197    echo "USAGE: sersic.norm (index)"
     1198    break
     1199  end
     1200
     1201  local index
     1202  $index = $1
     1203
     1204  if (($index >= 0.0) && ($index < 1.0))
     1205      $norm = 0.201545  - 0.950965 * $index - 0.315248 * $index^2
     1206      echo $norm {exp($norm)}
     1207      $myNorm = exp($norm)
     1208      return
     1209  end
     1210
     1211  if (($index >= 1.0) && ($index < 2.0))
     1212      $norm = 0.402084  - 1.357775 * $index - 0.105102 * $index^2
     1213      echo $norm {exp($norm)}
     1214      $myNorm = exp($norm)
     1215      return
     1216  end
     1217
     1218  if (($index >= 2.0) && ($index < 3.0))
     1219      $norm = 0.619093 - 1.591674 * $index - 0.041576 * $index^2
     1220      echo $norm {exp($norm)}
     1221      $myNorm = exp($norm)
     1222      return
     1223  end
     1224
     1225  if (($index >= 3.0) && ($index < 4.0))
     1226      $norm = 0.770263 - 1.696421 * $index - 0.023363 * $index^2
     1227      echo $norm {exp($norm)}
     1228      $myNorm = exp($norm)
     1229      return
     1230  end
     1231
     1232  if (($index >= 4.0) && ($index < 5.5))
     1233      $norm = 0.885891 - 1.755684 * $index - 0.015753 * $index^2
     1234      echo $norm {exp($norm)}
     1235      $myNorm = exp($norm)
     1236      return
     1237  end
     1238end
     1239
     1240macro load.model.im
     1241  if ($0 != 2)
     1242    echo "USAGE: load.model.im (N)"
     1243    break
     1244  end
     1245
     1246  rd obj$1 obj.$1.fits
     1247  rd cnv$1 cnv.$1.fits
     1248  rd var$1 var.$1.fits
     1249  rd msk$1 msk.$1.fits
     1250  for i 1 7
     1251    rd dpar$i.$1 dpar.$i.$1.fits
     1252  end
     1253  set dC$1 = (msk$1 == 0) * (obj$1 - cnv$1)^2 / var$1
     1254  tv dC$1 -0.01 3.0
     1255end
     1256
     1257macro load.normdata
     1258  if ($0 != 3)
     1259    echo "USAGE: load.normdata (file) [clear/noclear]"
    6641260    break
    6651261  end
    6661262
    6671263  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
     1264  read f 3 t 4 m 5
     1265  set M = -2.5*log(f)
     1266  subset Mg = M if (t == 1)
     1267  subset mg = m if (t == 1)
     1268  set dm = mg - Mg
     1269  set n = ramp(dm)
     1270  if ("$2" == "clear")
     1271    lim n dm; clear; box;
     1272  end
     1273  plot n dm
    6961274end
    6971275
    … …  
    7001278  exit 0
    7011279end
     1280
     1281# note that t1 = original confi
     1282# t2 = no residuals
     1283# t3 = const weight PSF fit
     1284# t4 = const weight galaxy fits
     1285
     1286macro run.radius.loop
     1287  local radius Nrun
     1288
     1289  if (1)
     1290    $Nrun = 0
     1291    foreach radius 1.0 1.5 2.0 2.5 3.0 4.0 6.0 8.0
     1292      mkexp.devexp.single tests.20131120/testrad.ps1v1.exp.$Nrun EXP $radius 1.0
     1293      fitexp tests.20131120/testrad.ps1v1.exp.$Nrun tests.20131120/testrad.ps1v1.exp.$Nrun.t4b EXP_CONV
     1294      $Nrun ++
     1295    end
     1296  end
     1297
     1298  $Nrun = 0
     1299  foreach radius 1.0 1.5 2.0 2.5 3.0 4.0 6.0 8.0
     1300    cmf.load.concat tests.20131120/testrad.ps1v1.exp.$Nrun.dat tests.20131120/testrad.ps1v1.exp.$Nrun.t4b.cmf EXP
     1301    echo -no-return $radius ""
     1302    check.fit tests.20131120/testrad.ps1v1.exp.$Nrun.dat tests.20131120/testrad.ps1v1.exp.$Nrun.t4b.cmf EXP
     1303    $Nrun ++
     1304  end
     1305
     1306  grid.plots.devexp a
     1307end
     1308
     1309macro check.radius.loop
     1310  if ($0 != 2)
     1311    echo "USAGE: check.radius.loop (ext)"
     1312    break
     1313  end
     1314
     1315  # cmf.load.reset
     1316
     1317  local radius Nrun
     1318
     1319  $Nrun = 0
     1320  foreach radius 1.0 1.5 2.0 2.5 3.0 4.0 6.0 8.0
     1321    echo -no-return $radius ""
     1322    check.fit tests.20131120/testrad.ps1v1.exp.$Nrun.dat tests.20131120/testrad.ps1v1.exp.$Nrun.$1.cmf EXP
     1323    $Nrun ++
     1324  end
     1325
     1326  grid.plots.devexp a
     1327end
     1328
     1329macro check.fit
     1330  if ($0 != 4)
     1331    echo "USAGE: check.fit (dat) (cmf) (type)"
     1332    break
     1333  end
     1334
     1335  cmf.load.concat $1 $2 $3
     1336  set dM_m = Mot_m - Min_m
     1337  vstat -q dM_m
     1338  $mag_off = $MEDIAN
     1339  set dR_m = Rot_m - Rin_m
     1340  vstat -q dR_m
     1341  $rad_off = $MEDIAN
     1342  echo $mag_off $rad_off
     1343end
     1344
     1345macro make.circles
     1346
     1347  delete Xo Yo Ro
     1348
     1349  for ix 0 16
     1350    for iy $ix 16
     1351      $r2 = $ix^2 + $iy^2
     1352      if ($r2 > 15^2) continue
     1353      concat $ix Xo
     1354      concat $iy Yo
     1355      concat $r2 Ro
     1356    end
     1357  end
     1358
     1359  sort Ro Xo Yo
     1360  $Rold = -1
     1361  $N = -1
     1362  for i 0 Ro[]
     1363    if (Ro[$i] != $Rold)
     1364      $N ++
     1365      $Rold = Ro[$i]
     1366    end
     1367    if ((Xo[$i] == 0) && (Yo[$i] == 0))
     1368      fprintf "// center is 0,0"
     1369      continue
     1370    end
     1371    if (Xo[$i] == 0)
     1372      fprintf "ADD_AXIS (%3d, %2d)     // r^2 = %3d" $N Yo[$i] Ro[$i]
     1373      continue
     1374    end
     1375    if (Xo[$i] == Yo[$i])
     1376      fprintf "ADD_DIAG (%3d, %2d)     // r^2 = %3d" $N Xo[$i] Ro[$i]
     1377      continue
     1378    end
     1379
     1380    fprintf "ADD_RAND (%3d, %2d, %2d) // r^2 = %3d" $N Xo[$i] Yo[$i] Ro[$i]
     1381 end
     1382end
     1383
     1384# EXP : tests.20131120/test.nsig
     1385# DEV : tests.20131120/test.dev
     1386
     1387# run.convolve.loop tests.20131120/test.nsig tests.20131120/test.nsig   EXP EXP_CONV
     1388# run.convolve.loop tests.20131120/test.dev  tests.20131120/test.dev    DEV DEV_CONV
     1389# run.convolve.loop tests.20131120/test.nsig tests.20131120/test.serexp EXP SER_CONV
     1390# run.convolve.loop tests.20131120/test.dev  tests.20131120/test.serdev DEV SER_CONV
     1391
     1392macro run.convolve.loop
     1393  if ($0 != 5)
     1394    echo "run.convolve.loop (inName) (fitName) (inType) (fitType)"
     1395    break
     1396  end
     1397
     1398  local radius
     1399  $radius = 3.0
     1400  cmf.load.reset
     1401
     1402  # foreach Cin 3 5 7 9 11
     1403  foreach Cin 11
     1404    $CONVOLVE_NSIGMA = $Cin
     1405    sprintf nameIn %s.%02d $1 $Cin
     1406    mkexp.devexp.single $nameIn $3 $radius 1.0
     1407    foreach Cot 3 5 7 9 11
     1408      $NSIGMA_CONV = $Cot
     1409      sprintf nameOt %s.%02d.%02d $2 $Cin $Cot
     1410      fitexp $nameIn $nameOt $4
     1411    end
     1412  end
     1413
     1414  # foreach Cin 3 5 7 9 11
     1415  foreach Cin 11
     1416    sprintf nameIn %s.%02d $1 $Cin
     1417    foreach Cot 3 5 7 9 11
     1418      sprintf nameOt %s.%02d.%02d $2 $Cin $Cot
     1419      check.fit $nameIn.dat $nameOt.cmf
     1420    end
     1421  end
     1422  grid.plots.devexp a
     1423end
     1424
     1425macro check.convolve.loop
     1426  if ($0 != 4)
     1427    echo "run.convolve.loop (inName) (fitName) (inType)"
     1428    break
     1429  end
     1430
     1431  local radius
     1432  $radius = 3.0
     1433  cmf.load.reset
     1434
     1435  # foreach Cin 3 5 7 9 11
     1436  foreach Cin 11
     1437    sprintf nameIn %s.%02d $1 $Cin
     1438    foreach Cot 3 5 7 9 11
     1439      sprintf nameOt %s.%02d.%02d $2 $Cin $Cot
     1440      echo -no-return $Cin $Cot " "
     1441      check.fit $nameIn.dat $nameOt.cmf $3
     1442    end
     1443  end
     1444  grid.plots.devexp a
     1445end
Note: See TracChangeset for help on using the changeset viewer.