IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 38280


Ignore:
Timestamp:
May 15, 2015, 1:56:42 PM (11 years ago)
Author:
eugene
Message:

fix interpolation of Sigma,Core -> FWHM in PS1_V1 and QGAUSS; test code to prove fix

Location:
trunk/psModules
Files:
3 added
6 edited

Legend:

Unmodified
Added
Removed

  • trunk/psModules/configure.ac

    r34085 r38280  
    77AC_CONFIG_SUBDIRS([test/tap])
    88
    9 AM_INIT_AUTOMAKE([1.7 foreign dist-bzip2])
     9AM_INIT_AUTOMAKE([1.7 foreign dist-bzip2 subdir-objects])
    1010AM_CONFIG_HEADER([src/config.h])
    1111AM_MAINTAINER_MODE

  • trunk/psModules/src/objects/models/fwhm.sh

    r36857 r38280  
    110110    $Zhm = 1.0
    111111    $FWHM = 2*sqrt(2)
    112     echo $K : $Zhm : $FWHM
     112    fprintf "%4.2f, // %4.1f, %4.2f" $FWHM $K $Zhm
    113113    return
    114114  end
    … …  
    121121 
    122122  set f = $K*z + z^1.667 - 1.0
    123   set dfdz = $K + z^0.667
     123  set dfdz = $K + 1.677*z^0.667
    124124
    125125  #lim -n 0 z f; clear; box; plot z f
    126126  #lim -n 1 z dfdz; clear; box; plot z dfdz 
    127127
    128   $nZ0 = 0
    129   $nZ1 = 5 / $dz
    130 
    131   $Zg = 0.5*(z[$nZ0] + z[$nZ1]) 
    132   $nZg = int($Zg / $dz)
    133   $dZ = 1.0
    134 
    135   for i 0 10
    136     $Fg = f[$nZg]
    137     $dFdz_g = dfdz[$nZg]
    138 
    139     $dZ = $Fg / $dFdz_g
    140 
    141     # echo $Zg $Fg $dZ $dFdz_g
    142 
    143     $Zg -= $dZ
    144     $Zg = max ($Zg , 0)
    145     $nZg = int($Zg / $dz)
     128  # these just need to bound the solution
     129  $nZ0 = 0.5 / $dz
     130  $nZ1 = 5.0 / $dz
     131
     132  $Zg = 0.5*(z[$nZ0] + z[$nZ1]) 
     133  $nZg = int($Zg / $dz)
     134  $dZ = 1.0
     135
     136  for i 0 10
     137    # interpolate between $nZg and $nZg + 1
     138    $Fg = f[$nZg] + ($Zg - $nZg*$dz)*(f[$nZg+1] - f[$nZg])/$dz
     139    $dFdz_g = dfdz[$nZg] + ($Zg - $nZg*$dz)*(dfdz[$nZg+1] - dfdz[$nZg])/$dz
     140
     141    # $Fg = f[$nZg]
     142    # $dFdz_g = dfdz[$nZg]
     143
     144    $dZ = $Fg / $dFdz_g
     145
     146    # echo -no-return $Zg $Fg $dZ $dFdz_g
     147
     148    $Zg -= $dZ
     149
     150    $Zg = max ($Zg , 0)
     151    $Zg = min ($Zg , z[-2])
     152    $nZg = int($Zg / $dz)
     153   
     154    # echo ": $Zg $nZg"
    146155  end
    147156  # echo $Zg $Fg $dZ $dFdz_g
    … …  
    149158  $FWHM = 2*sqrt(2*$Zg)
    150159
    151   echo $K : $Zhm : $FWHM
     160  fprintf "%4.2f, // %4.1f, %4.2f" $FWHM $K $Zhm
    152161end
    153162
    … …  
    164173    $Zhm = 1.0
    165174    $FWHM = 2*sqrt(2)
    166     echo $K : $Zhm : $FWHM
     175    fprintf "%4.2f, // %4.1f, %4.2f" $FWHM $K $Zhm
    167176    return
    168177  end
    … …  
    175184 
    176185  set f = $K*z + z^2.25 - 1.0
    177   set dfdz = $K + z^1.25
     186  set dfdz = $K + 2.25*z^1.25
    178187
    179188  #lim -n 0 z f; clear; box; plot z f
    … …  
    188197
    189198  for i 0 10
    190     $Fg = f[$nZg]
    191     $dFdz_g = dfdz[$nZg]
     199    # interpolate between $nZg and $nZg + 1
     200    $Fg = f[$nZg] + ($Zg - $nZg*$dz)*(f[$nZg+1] - f[$nZg])/$dz
     201    $dFdz_g = dfdz[$nZg] + ($Zg - $nZg*$dz)*(dfdz[$nZg+1] - dfdz[$nZg])/$dz
     202
     203    # $Fg = f[$nZg]
     204    # $dFdz_g = dfdz[$nZg]
    192205
    193206    $dZ = $Fg / $dFdz_g
    … …  
    197210    $Zg -= $dZ
    198211    $Zg = max ($Zg , 0)
     212    $Zg = min ($Zg , z[-2])
    199213    $nZg = int($Zg / $dz)
    200214  end
    … …  
    203217  $FWHM = 2*sqrt(2*$Zg)
    204218
    205   echo $K : $Zhm : $FWHM
     219  fprintf "%4.2f, // %4.1f, %4.2f" $FWHM $K $Zhm
    206220end
    207221
    208222macro fwhm.trend
     223  if ($0 != 4)
     224    echo "USAGE: fwhm.trend (model) (struct) (output)"
     225    echo "  model: pgauss, rgauss, ps1v1, qgauss"
     226    break
     227  end
    209228
    210229  delete fwhm_v k_v
    211230
    212   for k 0 20 1.0
    213     find.fwhm.qgauss $k
     231  $SAVE = 0
     232  if ($SAVE) exec rm -f $3
     233  if ($SAVE) output $3
     234
     235  $minK = -1
     236  $maxK = 20
     237  $delK = 0.2
     238
     239  if ($SAVE) echo "# define FWHM_BIN $delK"
     240  if ($SAVE) echo "# define MIN_FWHM_BIN $minK"
     241  if ($SAVE) echo "# define N_FWHM_BIN {int(($maxK - $minK) / $delK) + 1}"
     242  if ($SAVE) echo "static float $2[] = \{"
     243
     244  for k $minK $maxK $delK -incl
     245    find.fwhm.$1 $k
    214246    concat $k k_v
    215247    concat $FWHM fwhm_v
    216248  end
    217 end
    218 
     249  if ($SAVE) echo "\}@"
     250  if ($SAVE) output stdout
     251
     252  lim k_v fwhm_v; clear; box; plot k_v fwhm_v -pt 10 -sz 1.0
     253end
     254

  • trunk/psModules/src/objects/models/pmModel_PS1_V1.c

    r36857 r38280  
    337337}
    338338
    339 // I used the script in models/fwhm.sh to generate the trend of FWHM scaling vs the K value
    340 // FWHM = Scale * Sigma (not that PAR[PM_PAR_SXX] = sigma * sqrt(2)
    341 //  K : z_hm  : FWHM
    342 //  0 : 1.000 : 2.83
    343 //  1 : 0.597 : 2.19
    344 //  2 : 0.396 : 1.78
    345 //  3 : 0.291 : 1.53
    346 //  4 : 0.232 : 1.36
    347 //  5 : 0.198 : 1.26
    348 //  6 : 0.169 : 1.16
    349 //  7 : 0.142 : 1.07
    350 //  8 : 0.124 : 0.99
    351 //  9 : 0.118 : 0.97
    352 // 10 : 0.106 : 0.92
    353 // 11 : 0.092 : 0.86
    354 // 12 : 0.091 : 0.85
    355 // 13 : 0.080 : 0.80
    356 // 14 : 0.078 : 0.79
    357 // 15 : 0.073 : 0.76
    358 // 16 : 0.063 : 0.71
    359 // 17 : 0.068 : 0.74
    360 // 18 : 0.056 : 0.67
    361 // 19 : 0.058 : 0.68
    362 
    363 // static float PS1_V1_Core[] = { 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 17.0, 18.0, 19.0};
    364 static float PS1_V1_Scale[] = {2.83, 2.19, 1.78, 1.53, 1.36, 1.26, 1.16, 1.07, 0.99, 0.97, 0.92, 0.86, 0.85, 0.80, 0.79, 0.76, 0.71, 0.74, 0.67, 0.68};
     339// I used the script in models/fwhm.sh to generate the trend of FWHM scaling vs the K value:
     340# include "ps1v1.fwhm.h"
    365341
    366342psF64 PM_MODEL_SET_FWHM (const psVector *params, psF64 sigma) {
    … …  
    372348    if (!isfinite(core)) return (2.0*M_SQRT2*sigma);
    373349
    374     // if PS1_V1_Core is defined as a set of integer steps, so we can simplify:
    375     int binCore = MAX(0, MIN (19, (int)core));
    376 
    377     float scale = NAN;
    378     if (binCore == 0) {
    379         scale = (core - binCore + 0) * (PS1_V1_Scale[binCore + 1] - PS1_V1_Scale[binCore + 0]) + PS1_V1_Scale[binCore + 0];
    380     } else {
    381         scale = (core - binCore - 1) * (PS1_V1_Scale[binCore + 0] - PS1_V1_Scale[binCore - 1]) + PS1_V1_Scale[binCore - 1];
    382     }
     350    int binCore = MAX(0, MIN (N_FWHM_BIN - 2, (int)((core - MIN_FWHM_BIN)/FWHM_BIN)));
     351    float coreBin = binCore * FWHM_BIN + MIN_FWHM_BIN;
     352
     353    float scale = (core - coreBin) * (PS1_V1_Scale[binCore + 1] - PS1_V1_Scale[binCore]) / FWHM_BIN + PS1_V1_Scale[binCore];
     354
     355    // fprintf (stderr, "core: %f, binCore: %d, scale: %f\n", core, binCore, scale);
    383356
    384357    return (scale * sigma);

  • trunk/psModules/src/objects/models/pmModel_QGAUSS.c

    r36857 r38280  
    339339
    340340// I used the script in models/fwhm.sh to generate the trend of FWHM scaling vs the K value
    341 // FWHM = Scale * Sigma (not that PAR[PM_PAR_SXX] = sigma * sqrt(2)
    342 //  K : z_hm  : FWHM
    343 //  0 : 1.000 : 2.83
    344 //  1 : 0.648 : 2.28
    345 //  2 : 0.430 : 1.85
    346 //  3 : 0.310 : 1.58
    347 //  4 : 0.244 : 1.40
    348 //  5 : 0.200 : 1.26
    349 //  6 : 0.165 : 1.15
    350 //  7 : 0.149 : 1.09
    351 //  8 : 0.125 : 1.00
    352 //  9 : 0.116 : 0.96
    353 // 10 : 0.101 : 0.90
    354 // 11 : 0.095 : 0.87
    355 // 12 : 0.083 : 0.82
    356 // 13 : 0.082 : 0.81
    357 // 14 : 0.080 : 0.80
    358 // 15 : 0.074 : 0.77
    359 // 16 : 0.064 : 0.71
    360 // 17 : 0.068 : 0.74
    361 // 18 : 0.057 : 0.67
    362 // 19 : 0.058 : 0.68
    363 
    364 // static float QGAUSS_Core[] = { 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 17.0, 18.0, 19.0};
    365 static float QGAUSS_Scale[] = {2.83, 2.28, 1.85, 1.58, 1.40, 1.26, 1.15, 1.09, 1.00, 0.96, 0.90, 0.87, 0.82, 0.81, 0.80, 0.77, 0.71, 0.74, 0.67, 0.68};
     341# include "qgauss.fwhm.h"
    366342
    367343psF64 PM_MODEL_SET_FWHM (const psVector *params, psF64 sigma) {
    … …  
    373349    if (!isfinite(core)) return (2.0*M_SQRT2*sigma);
    374350
    375     // QGAUSS_Core is defined as a set of integer steps, so we can simplify:
    376     int binCore = MAX(0, MIN (19, (int)core));
    377 
    378     float scale = NAN;
    379     if (binCore == 0) {
    380         scale = (core - binCore + 0) * (QGAUSS_Scale[binCore + 1] - QGAUSS_Scale[binCore + 0]) + QGAUSS_Scale[binCore + 0];
    381     } else {
    382         scale = (core - binCore - 1) * (QGAUSS_Scale[binCore + 0] - QGAUSS_Scale[binCore - 1]) + QGAUSS_Scale[binCore - 1];
    383     }
     351    int binCore = MAX(0, MIN (N_FWHM_BIN - 2, (int)((core - MIN_FWHM_BIN)/FWHM_BIN)));
     352    float coreBin = binCore * FWHM_BIN + MIN_FWHM_BIN;
     353
     354    float scale = (core - coreBin) * (QGAUSS_Scale[binCore + 1] - QGAUSS_Scale[binCore]) / FWHM_BIN + QGAUSS_Scale[binCore];
    384355
    385356    return (scale * sigma);

  • trunk/psModules/test/objects/Makefile.am

    r36085 r38280  
    2121        tap_pmModel_CentralPixel \
    2222        tap_pmModel_CentralPixel_v2 \
     23        tap_pmModel_SET_FWHM \
    2324        tap_pmPSF \
    2425        tap_pmTrend2D \

  • trunk/psModules/test/objects/tap_pmPeaks.c

    r31153 r38280  
    218218                if (!((tmpPeak->type & PM_PEAK_LONE) || (tmpPeak->type & PM_PEAK_EDGE))) {
    219219                    diag("TEST ERROR: (0) peak at (%d, %d) (%f) ->type set improperly (0x%x).",
    220                           tmpPeak->y, tmpPeak->x, tmpPeak->value, tmpPeak->type);
     220                          tmpPeak->y, tmpPeak->x, tmpPeak->detValue, tmpPeak->type);
    221221                    diag(" should be (0x%x or 0x%x).\n", PM_PEAK_LONE, PM_PEAK_EDGE);
    222222                    testStatus = false;
    … …  
    225225                if (tmpPeak->type != PM_PEAK_LONE) {
    226226                    diag("TEST ERROR: (1) peak at (%d, %d) (%f) ->type set improperly (0x%x).\n",
    227                            tmpPeak->y, tmpPeak->x, tmpPeak->value, tmpPeak->type);
     227                           tmpPeak->y, tmpPeak->x, tmpPeak->detValue, tmpPeak->type);
    228228                    diag(" should be (0x%x).\n", PM_PEAK_LONE);
    229229                    testStatus = false;
    230230                }
    231231            } else {
    232                 diag("TEST ERROR: Peak at (%d, %d) (%f)\n", tmpPeak->y, tmpPeak->x, tmpPeak->value);
     232                diag("TEST ERROR: Peak at (%d, %d) (%f)\n", tmpPeak->y, tmpPeak->x, tmpPeak->detValue);
    233233                testStatus = false;
    234234            }
    … …  
    260260        ok(tmpPeak->x == 1, "pmPeakAlloc() set pmPeak->x");
    261261        ok(tmpPeak->y == 2, "pmPeakAlloc() set pmPeak->y");
    262         ok(tmpPeak->value == 3.0, "pmPeakAlloc() set pmPeak->value");
     262        ok(tmpPeak->detValue == 3.0, "pmPeakAlloc() set pmPeak->detValue");
    263263        ok(tmpPeak->flux == 0, "pmPeakAlloc() pmPeak->flux");
    264264        ok(tmpPeak->SN == 0, "pmPeakAlloc() pmPeak->SN");
Note: See TracChangeset for help on using the changeset viewer.