IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Dec 13, 2017, 10:53:48 AM (9 years ago)
Author:
eugene
Message:

merge EAM development branch changes for DR2 into trunk (add PS1_V6 dvo format; change Mcal to McalPSF, McalAPER; change opihi int vectors to 64bit)

Location:
trunk/Ohana
Files:
35 edited
2 copied

Legend:

Unmodified
Added
Removed

  • trunk/Ohana

  • trunk/Ohana/src/opihi/cmd.astro/fitplx.c

    r39610 r40291  
    7777  double *dD = dDvec->elements.Flt;
    7878
    79   int *mask = NULL;
     79  opihi_int *mask = NULL;
    8080  if (mvec) {
    8181    mask = mvec->elements.Int;
    … …  
    380380}
    381381
    382 int PlxSetMeanEpoch (double *R, double *D, double *T, double *Rmean, double *Dmean, double *Tmean, int *mask, int Ntotal) {
     382int PlxSetMeanEpoch (double *R, double *D, double *T, double *Rmean, double *Dmean, double *Tmean, opihi_int *mask, int Ntotal) {
    383383
    384384  int i;
    … …  
    413413
    414414// generate the fit values (projected X,Y; parallax factors;
    415 int PlxSetEpochPosition (PlxFitData *fitdata, double *R, double *D, double *dR, double *dD, double *T, int *mask, int Ntotal, Coords *coords, double Tmean) {
     415int PlxSetEpochPosition (PlxFitData *fitdata, double *R, double *D, double *dR, double *dD, double *T, opihi_int *mask, int Ntotal, Coords *coords, double Tmean) {
    416416
    417417  int i;
    … …  
    464464# define MAX_REJECT 0.1
    465465
    466 int PlxOutlierClip (PlxFitData *fitdata, int *mask, int Noutlier, float dPsigMax, Vector *dPvec, int VERBOSE) {
     466int PlxOutlierClip (PlxFitData *fitdata, opihi_int *mask, int Noutlier, float dPsigMax, Vector *dPvec, int VERBOSE) {
    467467
    468468  int i, n;

  • trunk/Ohana/src/opihi/cmd.astro/fitplx_irls.c

    r39926 r40291  
    8181  double *dD = dDvec->elements.Flt;
    8282
    83   int *mask = NULL;
     83  opihi_int *mask = NULL;
    8484  if (mvec) {
    8585    mask = mvec->elements.Int;
    … …  
    109109  for (i = 0; (VERBOSE == 2) && (i < fitdata.Npts); i++) {
    110110    int n = fitdata.index[i];
    111     int maskValue = mask ? mask[n] : 1;
    112     fprintf (stderr, "%f %f : %f %d : %f %f %f\n", R[n], D[n], T[n], maskValue, fitdata.t[i], fitdata.X[i], fitdata.Y[i]);
     111    opihi_int maskValue = mask ? mask[n] : 1;
     112    fprintf (stderr, "%f %f : %f "OPIHI_INT_FMT" : %f %f %f\n", R[n], D[n], T[n], maskValue, fitdata.t[i], fitdata.X[i], fitdata.Y[i]);
    113113  }
    114114
    … …  
    150150     
    151151      if (VERBOSE == 2) {
    152           fprintf (stderr, "%f %f : %f %d : %f %f %f : %f %f %f %f\n", R[n], D[n], T[n], mask[n], fitdata.t[i], fitdata.X[i], fitdata.Y[i], fitdata.Wx[i], fitdata.Wy[i], Sum_Wx, Sum_Wy);
     152          fprintf (stderr, "%f %f : %f "OPIHI_INT_FMT" : %f %f %f : %f %f %f %f\n", R[n], D[n], T[n], mask[n], fitdata.t[i], fitdata.X[i], fitdata.Y[i], fitdata.Wx[i], fitdata.Wy[i], Sum_Wx, Sum_Wy);
    153153      }
    154154    }

  • trunk/Ohana/src/opihi/cmd.astro/fitpm.c

    r39228 r40291  
    5050  double *dD = dDvec->elements.Flt;
    5151
    52   int *mask = NULL;
     52  opihi_int *mask = NULL;
    5353  if (mvec) {
    5454    mask = mvec->elements.Int;

  • trunk/Ohana/src/opihi/cmd.astro/fitpm_irls.c

    r39596 r40291  
    5858  double *dD = dDvec->elements.Flt;
    5959
    60   int *mask = NULL;
     60  opihi_int *mask = NULL;
    6161  if (mvec) {
    6262    mask = mvec->elements.Int;

  • trunk/Ohana/src/opihi/cmd.astro/star.c

    r36679 r40291  
    33int star (int argc, char **argv) {
    44
    5   int x, y, N, dx, Nborder;
     5  int x, y, N, Nborder;
    66  double max;
    77  Buffer *buf;
    … …  
    3333  }
    3434 
     35  int dx = 11;
     36  int dy = 11;
     37  int BOX = FALSE;
     38  if ((N = get_argument (argc, argv, "-box"))) {
     39    remove_argument (N, &argc, argv);
     40    dx  = atoi(argv[N]);
     41    remove_argument (N, &argc, argv);
     42    dy  = atoi(argv[N]);
     43    remove_argument (N, &argc, argv);
     44    BOX = TRUE;
     45  }
     46
    3547  if ((argc != 4) && (argc != 5)) {
    36     gprint (GP_ERR, "USAGE: star (buffer) x y [dx] [-border N] [-sat cnts]\n");
     48    gprint (GP_ERR, "USAGE: star (buffer) x y [dx] [-border N] [-sat cnts] [-box dx dy]\n");
    3749    gprint (GP_ERR, " dx is the aperture diameter, but is adjusted up to the next odd number\n");
    3850    return (FALSE);
    … …  
    4052  if ((buf = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) return (FALSE);
    4153
    42   dx = 11;
    4354  x = atof (argv[2]);
    4455  y = atof (argv[3]);
    … …  
    4758  }
    4859
    49   get_aperture_stats (&buf[0].matrix, x, y, dx, Nborder, max, VERBOSE);
     60  if (BOX) {
     61    get_box_stats (&buf[0].matrix, x, y, dx, dy, Nborder, max, VERBOSE);
     62  } else {
     63    get_aperture_stats (&buf[0].matrix, x, y, dx, Nborder, max, VERBOSE);
     64  }
    5065 
    5166  return (TRUE);

  • trunk/Ohana/src/opihi/cmd.data/Makefile

    r40165 r40291  
    155155$(SRC)/type.$(ARCH).o              \
    156156$(SRC)/uniq.$(ARCH).o              \
     157$(SRC)/uniqpair.$(ARCH).o                  \
    157158$(SRC)/unsign.$(ARCH).o            \
    158159$(SRC)/vbin.$(ARCH).o              \

  • trunk/Ohana/src/opihi/cmd.data/init.c

    r40165 r40291  
    141141int tvcontour        PROTO((int, char **));
    142142int tvgrid           PROTO((int, char **));
    143 int opihi_type             PROTO((int, char **));
     143int opihi_type       PROTO((int, char **));
    144144int uniq             PROTO((int, char **));
     145int uniqpair         PROTO((int, char **));
    145146int unsign           PROTO((int, char **));
    146147int vbin             PROTO((int, char **));
    … …  
    325326  {1, "ungridify",    ungridify,        "convert image region to vector triplet"},
    326327  {1, "uniq",         uniq,             "create a uniq vector subset from a vector"},
     328  {1, "uniqpair",     uniqpair,         "create a uniq vector subset from a pair of vectors, saving duplicates if desired"},
    327329  {1, "unsign",       unsign,           "toggle the UNSIGN status"},
    328330  {1, "vbin",         vbin,             "rebin vector data by a factor of N"},

  • trunk/Ohana/src/opihi/cmd.data/limits.c

    r31160 r40291  
    33int limits (int argc, char **argv) {
    44
    5   int N, APPLY, dX, dY;
     5  int N, dX, dY;
    66  int kapa;
    7   char *name;
    87  Graphdata graphmode;
    98  Vector *xvec, *yvec;
    … …  
    1110  xvec = yvec = NULL;
    1211
    13   APPLY = FALSE;
     12  float minLimitX = NAN;
     13  float minLimitY = NAN;
     14  float maxLimitX = NAN;
     15  float maxLimitY = NAN;
     16  float delLimitX = NAN;
     17  float delLimitY = NAN;
     18
     19  if ((N = get_argument (argc, argv, "-minX"))) {
     20    remove_argument (N, &argc, argv);
     21    minLimitX = atof (argv[N]);
     22    remove_argument (N, &argc, argv);
     23  }
     24  if ((N = get_argument (argc, argv, "-maxX"))) {
     25    remove_argument (N, &argc, argv);
     26    maxLimitX = atof (argv[N]);
     27    remove_argument (N, &argc, argv);
     28  }
     29  if ((N = get_argument (argc, argv, "-delX"))) {
     30    if (!isnan(minLimitX) || !isnan(maxLimitX)) {
     31      gprint (GP_ERR, "-minX & -maxX cannot be mixed with -delX\n");
     32      return (FALSE);
     33    }
     34    remove_argument (N, &argc, argv);
     35    delLimitX = atof (argv[N]);
     36    remove_argument (N, &argc, argv);
     37  }
     38  if ((N = get_argument (argc, argv, "-minY"))) {
     39    remove_argument (N, &argc, argv);
     40    minLimitY = atof (argv[N]);
     41    remove_argument (N, &argc, argv);
     42  }
     43  if ((N = get_argument (argc, argv, "-maxY"))) {
     44    remove_argument (N, &argc, argv);
     45    maxLimitY = atof (argv[N]);
     46    remove_argument (N, &argc, argv);
     47  }
     48  if ((N = get_argument (argc, argv, "-delY"))) {
     49    if (!isnan(minLimitY) || !isnan(maxLimitY)) {
     50      gprint (GP_ERR, "-minY & -maxY cannot be mixed with -delY\n");
     51      return (FALSE);
     52    }
     53    remove_argument (N, &argc, argv);
     54    delLimitY = atof (argv[N]);
     55    remove_argument (N, &argc, argv);
     56  }
     57
     58  int APPLY = FALSE;
    1459  if ((N = get_argument (argc, argv, "-a"))) {
    1560    remove_argument (N, &argc, argv);
    1661    APPLY = TRUE;
    1762  }
    18   name = NULL;
     63  char *name = NULL;
    1964  if ((N = get_argument (argc, argv, "-n"))) {
    2065    remove_argument (N, &argc, argv);
    … …  
    2267    remove_argument (N, &argc, argv);
    2368  }
     69
    2470  if (!GetGraph (&graphmode, &kapa, name)) return (FALSE);
    2571  FREE (name);
    … …  
    97143 success:
    98144  SetLimits (xvec, yvec, &graphmode);
     145
     146  if (!isnan(minLimitX)) graphmode.xmin = MIN (minLimitX, graphmode.xmin);
     147  if (!isnan(maxLimitX)) graphmode.xmax = MAX (maxLimitX, graphmode.xmax);
     148  if (!isnan(minLimitY)) graphmode.ymin = MIN (minLimitY, graphmode.ymin);
     149  if (!isnan(maxLimitY)) graphmode.ymax = MAX (maxLimitY, graphmode.ymax);
     150
     151  if (!isnan(delLimitX)) {
     152    float delta = graphmode.xmax - graphmode.xmin;
     153    if (fabs(delLimitX) > fabs(delta)) {
     154      float midpt = 0.5*(graphmode.xmax + graphmode.xmin);
     155      graphmode.xmax = midpt + 0.5*delLimitX;
     156      graphmode.xmin = midpt - 0.5*delLimitX;
     157    }
     158  }
     159  if (!isnan(delLimitY)) {
     160    float delta = graphmode.ymax - graphmode.ymin;
     161    if (fabs(delLimitY) > fabs(delta)) {
     162      float midpt = 0.5*(graphmode.ymax + graphmode.ymin);
     163      graphmode.ymax = midpt + 0.5*delLimitY;
     164      graphmode.ymin = midpt - 0.5*delLimitY;
     165    }
     166  }
     167
    99168  if (APPLY) KapaSetLimits (kapa, &graphmode);
    100169  return (TRUE);

  • trunk/Ohana/src/opihi/cmd.data/print_vectors.c

    r37049 r40291  
    44
    55  Vector **vec;
    6   int i, j;
     6  int i, j, N;
     7
     8  int START_VALUE = 0;
     9  if ((N = get_argument (argc, argv, "-s"))) {
     10    remove_argument (N, &argc, argv);
     11    START_VALUE = atoi (argv[N]);
     12    remove_argument (N, &argc, argv);
     13  }
     14
     15  int END_VALUE = -1;
     16  if ((N = get_argument (argc, argv, "-e"))) {
     17    remove_argument (N, &argc, argv);
     18    END_VALUE = atoi (argv[N]);
     19    remove_argument (N, &argc, argv);
     20  }
    721
    822  if (argc < 2) {
    … …  
    2741  }
    2842
    29   for (j = 0; j < MaxLen; j++) {
     43  // start and end may be 0 - N (truncated to N) or may be negative, in which case it refers to
     44  // distance from the end (just like vector[-5])
     45  START_VALUE = (START_VALUE < 0) ? MaxLen + START_VALUE + 1 : MIN (START_VALUE, MaxLen);
     46  START_VALUE = MAX (0, START_VALUE);
     47
     48  END_VALUE = (END_VALUE < 0) ? MaxLen + END_VALUE + 1 : MIN (END_VALUE, MaxLen);
     49  END_VALUE = MAX (0, END_VALUE);
     50
     51  for (j = START_VALUE; j < END_VALUE; j++) {
    3052    for (i = 0; i < Nvec; i++) {
    3153      if (j >= vec[i][0].Nelements) {
    … …  
    3557          gprint (GP_LOG, "%f ", vec[i][0].elements.Flt[j]);
    3658        } else {
    37           gprint (GP_LOG, "%d ", vec[i][0].elements.Int[j]);
     59          gprint (GP_LOG, OPIHI_INT_FMT" ", vec[i][0].elements.Int[j]);
    3860        }
    3961      }

  • trunk/Ohana/src/opihi/cmd.data/reindex.c

    r39227 r40291  
    4747        continue;
    4848      }
    49       if (*vx > Nmax) ESCAPE("unexpected value in index: %d (%d)\n", *vx, i);
     49      if (*vx > Nmax) ESCAPE("unexpected value in index: "OPIHI_INT_FMT" (%d)\n", *vx, i);
    5050      ovec[0].elements.Flt[Npts] = vi[*vx];
    5151      Npts++;
    … …  
    6767        continue;
    6868      }
    69       if (*vx > Nmax) ESCAPE("unexpected value in index: %d (%d)\n", *vx, i);
     69      if (*vx > Nmax) ESCAPE("unexpected value in index: "OPIHI_INT_FMT" (%d)\n", *vx, i);
    7070      ovec[0].elements.Int[Npts] = vi[*vx];
    7171      Npts++;

  • trunk/Ohana/src/opihi/cmd.data/test/periodogram-fm.sh

    r40165 r40291  
    5555
    5656 periodogram_fm t f df 5 50 period power
    57 #periodogram t f 5 50 period power
    5857
    5958 # lim -n 0 t f; clear; box; plot -x 2 -pt 2 t f
    … …  
    8887
    8988 periodogram_fm t f df 1 10 period power
    90 #periodogram t f 1 10 period power
    9189
    9290 # lim -n 0 t f; clear; box; plot -x 2 -pt 2 t f
    … …  
    9795 if (abs ($peakpos - $P) > 0.05)
    9896   $PASS = 0
     97 end
     98
     99  if ($PLOT)
     100  lim period power; clear; box; line -c red70 -lw 3 $P 0 to $P $peakval; plot period power -x line
    99101 end
    100102end
    … …  
    117119
    118120 periodogram_fm t f df 2 30 period power
    119 #periodogram t f 2 30 period power
    120121
    121122#  lim -n 0 t f; clear; box; plot -x 2 -pt 2 t f
    … …  
    126127 if (abs ($peakpos - $P) > 0.05)
    127128   $PASS = 0
     129 end
     130
     131 if ($PLOT)
     132  lim period power; clear; box; line -c red70 -lw 3 $P 0 to $P $peakval; plot period power -x line
    128133 end
    129134end
    … …  
    145150 set df = 0.01 + zero(f)
    146151
    147  periodogram_fm t f 2 30 period power
     152 periodogram_fm t f df 2 30 period power
    148153
    149154#  lim -n 0 t f; clear; box; plot -x 2 -pt 2 t f
    … …  
    155160   $PASS = 0
    156161 end
    157 end
    158 
    159 # test using random samples, offset start, non-zero DC, some noise
     162
     163 if ($PLOT)
     164  lim period power; clear; box; line -c red70 -lw 3 $P 0 to $P $peakval; plot period power -x line
     165 end
     166end
     167
     168# test using 300 random samples, offset start, non-zero DC, some noise
    160169macro test6
    161170 $PASS = 1
    … …  
    178187 set df = 0.01 + zero(f)
    179188
    180  periodogram_fm t f 2 30 period power
     189 periodogram_fm t f df 2 30 period power
    181190
    182191#  lim -n 0 t f; clear; box; plot -x 2 -pt 2 t f
    … …  
    188197   $PASS = 0
    189198 end
    190 end
    191 
    192 # test using fewer random samples, offset start, non-zero DC, some noise
     199
     200 if ($PLOT)
     201  lim period power; clear; box; line -c red70 -lw 3 $P 0 to $P $peakval; plot period power -x line
     202 end
     203end
     204
     205# test using 100 fewer random samples, offset start, non-zero DC, some noise
    193206macro test7
    194207 $PASS = 1
    … …  
    211224 set df = 0.01 + zero(f)
    212225
    213  periodogram_fm t f 2 30 period power
     226 periodogram_fm t f df 2 30 period power
    214227
    215228#  lim -n 0 t f; clear; box; plot -x 2 -pt 2 t f
    … …  
    221234   $PASS = 0
    222235 end
    223 end
    224 
    225 # test using fewer random samples, high frequency, non-zero DC, some noise
     236
     237 if ($PLOT)
     238  lim period power; clear; box; line -c red70 -lw 3 $P 0 to $P $peakval; plot period power -x line
     239 end
     240end
     241
     242# test using Ndays random samples, RR Lyrae-sized light curves (0.7 mag),
     243# optional noise level
    226244macro test8
    227  if ($0 != 2)
    228    echo "USAGE: test8: Ndays")
     245 if ($0 != 4)
     246   echo "USAGE: test8: Period Ndays (df)"
    229247   break
    230248 end
    231249 
    232250 local Ndays
    233  $Ndays = $1
    234 
    235  $PASS = 1
    236  break -auto off
    237 
    238  local P PI
    239  $PI = 3.14159265359
    240  $P  = 0.8*rnd(0) + 0.2
     251 $P = $1
     252 $Ndays = $2
     253 $dM = $3
     254
     255 $PASS = 1
     256 break -auto off
     257
     258 local PI
     259 $PI = 3.14159265359
    241260 $trueP = $P
    242261
    … …  
    246265
    247266 # t is a time in days, but we always have 4 within 1 hour:
    248  set t0 = int(100 * rnd(x))
    249  set dtx = (3/24) * rnd(x)
    250  set t0 = t0 + dtx
     267 set tday = int(100 * rnd(x)); # choose Ndays random days between 0 and 100
     268 set dtx = (3/24) * rnd(x);  # choose a starting time within that night
     269 set t0 = tday + dtx
    251270
    252271 set dt1 = (15.0 / 1440) * rnd(x) + ( 0 + 7.5) / 1440
    … …  
    260279 set tmp = t0 + dt3; concat tmp t
    261280
    262  set fraw = sin(2*$PI*t/$P) + 0.5
     281 set fraw = 0.75*sin(2*$PI*t/$P)
    263282
    264283 # 0.05 : peakpos = 14.95
    265284 # 0.10 : peakpos = 15.04 (
    266  gaussdev df t[] 0.0 0.25
     285 gaussdev df t[] 0.0 $dM
    267286 set f = fraw + df
    268  set df = 0.01 + zero(f)
    269 
    270  periodogram_fm t f 0.1 2.0 period power
     287 set df = $dM + zero(f)
     288
     289 periodogram_fm t f df 0.1 20.0 period power
    271290
    272291#  lim -n 0 t f; clear; box; plot -x 2 -pt 2 t f
    … …  
    278297   $PASS = 0
    279298 end
    280 end
     299 if ($PLOT)
     300  lim period power; clear; box; line -c red70 -lw 3 $P 0 to $P $peakval; plot period power -x line
     301 end
     302end
     303
     304# test using Ndays random samples, RR Lyrae-sized light curves (0.7 mag),
     305# optional noise level
     306# compare periodogram and periodogram_fm
     307macro test9
     308 if ($0 != 4)
     309   echo "USAGE: test8: Period Ndays (df)"
     310   break
     311 end
     312 
     313 local Ndays
     314 $P = $1
     315 $Ndays = $2
     316 $dM = $3
     317
     318 $PASS = 1
     319 break -auto off
     320
     321 local PI
     322 $PI = 3.14159265359
     323 $trueP = $P
     324
     325 delete -q x t f period power
     326
     327 create x 0 $Ndays
     328
     329 # t is a time in days, but we always have 4 within 1 hour:
     330 set tday = int(100 * rnd(x)); # choose Ndays random days between 0 and 100
     331 set dtx = (3/24) * rnd(x);  # choose a starting time within that night
     332 set t0 = tday + dtx
     333
     334 set dt1 = (15.0 / 1440) * rnd(x) + ( 0 + 7.5) / 1440
     335 set dt2 = (15.0 / 1440) * rnd(x) + (15 + 7.5) / 1440
     336 set dt3 = (15.0 / 1440) * rnd(x) + (30 + 7.5) / 1440
     337
     338 delete -q t
     339 concat t0 t
     340 set tmp = t0 + dt1; concat tmp t
     341 set tmp = t0 + dt2; concat tmp t
     342 set tmp = t0 + dt3; concat tmp t
     343
     344 set fraw = 0.75*sin(2*$PI*t/$P)
     345
     346 # 0.05 : peakpos = 14.95
     347 # 0.10 : peakpos = 15.04 (
     348 gaussdev df t[] 0.0 $dM
     349 set f = fraw + df
     350 set df = $dM + zero(f)
     351
     352 periodogram_fm t f df 0.1 20.0 period_fm power_fm
     353 periodogram t f 0.1 20.0 period power
     354
     355#  lim -n 0 t f; clear; box; plot -x 2 -pt 2 t f
     356#  lim -n 1 period power; clear; box; plot period power
     357
     358 peak -q period_fm power_fm
     359 $peakval_fm = $peakval
     360
     361 peak -q period power
     362 vstat -q power
     363 set power = power / $MAX
     364
     365# if (abs ($peakpos - $P) > 0.05)
     366#   $PASS = 0
     367# end
     368
     369 set freq = 1 / period
     370 set freq_fm = 1 / period_fm
     371 $Freq = 1 / $P
     372
     373 if ($PLOT)
     374  if (1)
     375    lim period power; clear; box
     376    line -c red70 -lw 3 $P 0 to $P $peakval_fm;
     377    plot period power -x line -c grey70 -lw 2
     378    plot period_fm power_fm -x line -c black
     379  else
     380    lim freq power; clear; box
     381    line -c red70 -lw 3 $Freq 0 to $Freq 1.0
     382    plot freq power -x line -c grey70 -lw 2
     383    plot freq_fm power_fm -x line -c black
     384  end
     385 end
     386end
     387
     388
     389# test using Ndays random samples, RR Lyrae-sized light curves (0.7 mag),
     390# optional noise level
     391# compare periodogram and periodogram_fm
     392macro test10
     393 if ($0 != 4)
     394   echo "USAGE: test8: Period Ndays (df)"
     395   break
     396 end
     397 
     398 local Ndays
     399 $P = $1
     400 $Ndays = $2
     401 $dM = $3
     402
     403 $PASS = 1
     404 break -auto off
     405
     406 local PI
     407 $PI = 3.14159265359
     408 $trueP = $P
     409
     410 delete -q x t f period power
     411
     412 create x 0 $Ndays
     413
     414 # t is a time in days, but we always have 4 within 1 hour:
     415 set tday = int(100 * rnd(x)); # choose Ndays random days between 0 and 100
     416 set dtx = (3/24) * rnd(x);  # choose a starting time within that night
     417 set t0 = tday + dtx
     418
     419 set dt1 = (15.0 / 1440) * rnd(x) + ( 0 + 7.5) / 1440
     420 set dt2 = (15.0 / 1440) * rnd(x) + (15 + 7.5) / 1440
     421 set dt3 = (15.0 / 1440) * rnd(x) + (30 + 7.5) / 1440
     422
     423 delete -q t
     424 concat t0 t
     425 set tmp = t0 + dt1; concat tmp t
     426 set tmp = t0 + dt2; concat tmp t
     427 set tmp = t0 + dt3; concat tmp t
     428
     429 set fraw = 0.75*sin(2*$PI*t/$P)
     430
     431 # 0.05 : peakpos = 14.95
     432 # 0.10 : peakpos = 15.04 (
     433 gaussdev df t[] 0.0 $dM
     434 set f = fraw + df
     435 set df = $dM + zero(f)
     436
     437 periodogram_fm t f df 0.05 20.0 period_fm power_fm
     438
     439 gaussdev df t[] 0.0 $dM
     440 set Fo = df
     441 periodogram_fm t Fo df 0.05 20.0 period power
     442
     443#  lim -n 0 t f; clear; box; plot -x 2 -pt 2 t f
     444#  lim -n 1 period power; clear; box; plot period power
     445
     446 peak -q period_fm power_fm
     447 $peakval_fm = $peakval
     448
     449 peak -q period power
     450
     451# if (abs ($peakpos - $P) > 0.05)
     452#   $PASS = 0
     453# end
     454
     455 set freq = 1 / period
     456 set freq_fm = 1 / period_fm
     457 $Freq = 1 / $P
     458
     459 if ($PLOT)
     460  if (1)
     461    lim period_fm power_fm; clear; box
     462    line -c red70 -lw 3 $P 0 to $P $peakval_fm;
     463    plot period power -x line -c grey70 -lw 2
     464    plot period_fm power_fm -x line -c black
     465  else
     466    lim freq power; clear; box
     467    line -c red70 -lw 3 $Freq 0 to $Freq 1.0
     468    plot freq power -x line -c grey70 -lw 2
     469    plot freq_fm power_fm -x line -c black
     470  end
     471 end
     472end
     473
     474# we have time (MJD) and mag
     475# we generate the folded lightcure and measure sigma relative to the smoothed version (bins of 0.1 period)
     476macro fold.one.period
     477  if ($0 != 5)
     478    echo "USAGE: fold.one.period (time) (mag) (magErr) (period)"
     479    break
     480  end
     481
     482  local myTime myMag myMagErr myPeriod
     483  $myTime = $1
     484  $myMag  = $2
     485  $myMagErr  = $3
     486  $myPeriod = $4
     487
     488  set phi = $myTime / $myPeriod - int($myTime / $myPeriod)
     489
     490  if ($PLOT_FOLD)
     491    lim -n phi phi $myMag; clear; box;
     492  end
     493
     494  delete -q magResid
     495
     496  $dPhi = 0.05; # half of bin size
     497  create nphi $dPhi {1 + $dPhi} {2*$dPhi}
     498  set magR = zero(nphi)
     499  set magS = zero(nphi)
     500  for i 0 nphi[]
     501    subset tmp_mag_sub = $myMag where (phi >= nphi[$i] - $dPhi) && (phi < nphi[$i] + $dPhi)
     502    vstat -q tmp_mag_sub
     503    magR[$i] = $MEDIAN
     504    magS[$i] = $SIGMA
     505
     506    set magDelta = tmp_mag_sub - $MEDIAN
     507    concat magDelta magResid
     508
     509    if ($PLOT_FOLD)
     510      subset tmp_phi_sub = phi where (phi >= nphi[$i] - $dPhi) && (phi < nphi[$i] + $dPhi)
     511      if ($i % 2)
     512        plot tmp_phi_sub tmp_mag_sub -pt 7 -sz 3 -c blue -lw 2
     513      else
     514        plot tmp_phi_sub tmp_mag_sub -pt 7 -sz 3 -c red -lw 2
     515      end
     516    end 
     517  end
     518
     519  if ($PLOT_FOLD) 
     520    plot -pt 10 -sz 1.5 phi $myMag -dy $myMagErr
     521    plot -pt 2 -sz 2.0 -c red nphi magR -dy magS
     522  end
     523
     524  vstat -q magResid
     525end
     526
     527
    281528
    282529# Memory test
    … …  
    298545
    299546 for i 0 100
    300   periodogram_fm t f 2 30 period power
     547  periodogram_fm t f df 2 30 period power
    301548 end
    302549 

  • trunk/Ohana/src/opihi/cmd.data/test/periodogram.sh

    r40165 r40291  
    209209# test using fewer random samples, high frequency, non-zero DC, some noise
    210210macro test8
    211  if ($0 != 2)
    212    echo "USAGE: test8: Ndays")
     211 if ($0 != 4)
     212   echo "USAGE: test8: Period Ndays df"
    213213   break
    214214 end
    215215 
    216216 local Ndays
    217  $Ndays = $1
    218 
    219  $PASS = 1
    220  break -auto off
    221 
    222  local P PI
    223  $PI = 3.14159265359
    224  $P  = 0.8*rnd(0) + 0.2
     217 $P = $1
     218 $Ndays = $2
     219 $dM = $3
     220
     221 $PASS = 1
     222 break -auto off
     223
     224 local PI
     225 $PI = 3.14159265359
    225226 $trueP = $P
    226227
    … …  
    230231
    231232 # t is a time in days, but we always have 4 within 1 hour:
    232  set t0 = int(100 * rnd(x))
    233  set dtx = (3/24) * rnd(x)
    234  set t0 = t0 + dtx
     233 set tday = int(100 * rnd(x)); # choose Ndays random days between 0 and 100
     234 set dtx = (3/24) * rnd(x);  # choose a starting time within that night
     235 set t0 = tday + dtx
    235236
    236237 set dt1 = (15.0 / 1440) * rnd(x) + ( 0 + 7.5) / 1440
    … …  
    244245 set tmp = t0 + dt3; concat tmp t
    245246
    246  set fraw = sin(2*$PI*t/$P) + 0.5
     247 set fraw = 0.75*sin(2*$PI*t/$P)
    247248
    248249 # 0.05 : peakpos = 14.95
    249250 # 0.10 : peakpos = 15.04 (
    250  gaussdev df t[] 0.0 0.25
     251 gaussdev df t[] 0.0 $dM
    251252 set f = fraw + df
    252253
    … …  
    260261 if (abs ($peakpos - $P) > 0.05)
    261262   $PASS = 0
     263 end
     264 if ($PLOT)
     265  lim period power; clear; box; line -c red70 -lw 3 $P 0 to $P $peakval; plot period power -x line
    262266 end
    263267end

  • trunk/Ohana/src/opihi/cmd.data/uniq.c

    r39457 r40291  
    8080    memcpy (indata, ivec->elements.Int, ivec[0].Nelements*sizeof(opihi_int));
    8181
    82     isort (indata, ivec->Nelements);
     82    llsort (indata, ivec->Nelements);
    8383
    8484    Nnew = 0;

  • trunk/Ohana/src/opihi/cmd.data/write_vectors.c

    r39360 r40291  
    165165        } else {
    166166          if (CSV) {
    167             fprintf (f, "%d,", vec[j][0].elements.Int[i]);
     167            fprintf (f, OPIHI_INT_FMT",", vec[j][0].elements.Int[i]);
    168168          } else {
    169             fprintf (f, "%d ", vec[j][0].elements.Int[i]);
     169            fprintf (f, OPIHI_INT_FMT" ", vec[j][0].elements.Int[i]);
    170170          }
    171171        }

  • trunk/Ohana/src/opihi/dvo/dvo_host_utils.c

    r39524 r40291  
    284284    // XXX a bit of a waste (but only 1024 * 60 bytes or so
    285285    ALLOCATE (table->hosts[i].results, char, DVO_MAX_PATH);
    286     snprintf (table->hosts[i].results, DVO_MAX_PATH, "%s/dvo.results.%s.fits", table->hosts[i].pathname, uniquer);
     286    snprintf (table->hosts[i].results, DVO_MAX_PATH, "%s/dvo.results.%s.%04d.fits", table->hosts[i].pathname, uniquer, table->hosts[i].hostID);
    287287
    288288    int    Ninvec = 0;

  • trunk/Ohana/src/opihi/dvo/gimages.c

    r39347 r40291  
    214214    if (PixelCoords) {
    215215      gprint (GP_LOG, "%3d %5d %s %6.1f %6.1f %20s %5d %2d %4.2f %6.3f %5.3f %5.3f %4x %7d\n",
    216               Nfound, (int) i, image[i].name, X, Y, date, image[i].nstar, image[i].photcode, image[i].secz, image[i].Mcal, image[i].dMcal, image[i].exptime, image[i].flags, image[i].imageID);
     216              Nfound, (int) i, image[i].name, X, Y, date, image[i].nstar, image[i].photcode, image[i].secz, image[i].McalPSF, image[i].dMcal, image[i].exptime, image[i].flags, image[i].imageID);
    217217    } else {
    218218      XY_to_RD (&ra, &dec, 0.5*image[i].NX, 0.5*image[i].NY, &image[i].coords);
    219219      gprint (GP_LOG, "%3d %5d %s %8.4f %8.4f %20s %5d %2d %4.2f %6.3f %5.3f %5.3f %4x %7d\n",
    220               Nfound, (int) i, image[i].name, ra, dec, date, image[i].nstar, image[i].photcode, image[i].secz, image[i].Mcal, image[i].dMcal, image[i].exptime, image[i].flags, image[i].imageID);
     220              Nfound, (int) i, image[i].name, ra, dec, date, image[i].nstar, image[i].photcode, image[i].secz, image[i].McalPSF, image[i].dMcal, image[i].exptime, image[i].flags, image[i].imageID);
    221221    }
    222222    sprintf (name, "IMAGEx:%d", Nfound);

  • trunk/Ohana/src/opihi/dvo/gstar.c

    r39634 r40291  
    760760
    761761            if (FULL_OUTPUT) {
    762               gprint (GP_LOG, "%6.3f ", catalog.measure[Nv].Mcal);
     762              gprint (GP_LOG, "%6.3f ", catalog.measure[Nv].McalPSF);
     763              gprint (GP_LOG, "%6.3f ", catalog.measure[Nv].McalAPER);
    763764              gprint (GP_LOG, "%6.3f ", catalog.measure[Nv].Mflat);
    764               gprint (GP_LOG, "%6.3f ", catalog.measure[Nv].Map);
    765               gprint (GP_LOG, "%6.3f ", catalog.measure[Nv].Mkron);
     765              Mrel = PhotRel (&catalog.measure[Nv], &catalog.average[k], &catalog.secfilt[k*Nsecfilt], MAG_CLASS_APER);
     766              gprint (GP_LOG, "%6.3f ", Mrel);
     767              Mrel = PhotRel (&catalog.measure[Nv], &catalog.average[k], &catalog.secfilt[k*Nsecfilt], MAG_CLASS_KRON);
     768              gprint (GP_LOG, "%6.3f ", Mrel);
    766769              gprint (GP_LOG, "%6.3f ", catalog.measure[Nv].dMkron);
    767770              gprint (GP_LOG, "%5.1f ", pow(10.0, 0.4*catalog.measure[Nv].dt));
    … …  
    960963        print_double (NAN);
    961964      } else {
    962         print_double_exp (secfilt[seq].Mstdev);
     965        print_double_exp (secfilt[seq].sMpsfChp);
    963966      }
    964967      break;
    … …  
    10231026        print_double (NAN);
    10241027      } else {
    1025         print_double (secfilt[seq].M);
     1028        print_double (secfilt[seq].MpsfChp);
    10261029      }
    10271030      break;
    … …  
    10311034        print_double (NAN);
    10321035      } else {
    1033         print_double (secfilt[seq].dM);
     1036        print_double (secfilt[seq].dMpsfChp);
    10341037      }
    10351038      break;
    … …  
    10391042        print_double (NAN);
    10401043      } else {
    1041         print_double (secfilt[seq].Map);
     1044        print_double (secfilt[seq].MapChp);
    10421045      }
    10431046      break;
    … …  
    10471050        print_double (NAN);
    10481051      } else {
    1049         print_double (secfilt[seq].dMap);
     1052        print_double (secfilt[seq].dMapChp);
    10501053      }
    10511054      break;
    … …  
    10551058        print_double (NAN);
    10561059      } else {
    1057         print_double_exp (secfilt[seq].sMap);
     1060        print_double_exp (secfilt[seq].sMapChp);
    10581061      }
    10591062      break;
    … …  
    10631066        print_double (NAN);
    10641067      } else {
    1065         print_double_exp (secfilt[seq].dMap);
     1068        print_double_exp (secfilt[seq].dMapChp);
    10661069      }
    10671070      break;
    … …  
    10711074        print_double (NAN);
    10721075      } else {
    1073         print_double (secfilt[seq].Mkron);
     1076        print_double (secfilt[seq].MkronChp);
    10741077      }
    10751078      break;
    … …  
    10791082        print_double (NAN);
    10801083      } else {
    1081         print_double (secfilt[seq].dMkron);
     1084        print_double (secfilt[seq].dMkronChp);
    10821085      }
    10831086      break;
    … …  
    10871090        print_double (NAN);
    10881091      } else {
    1089         print_double_exp (secfilt[seq].sMkron);
     1092        print_double_exp (secfilt[seq].sMkronChp);
    10901093      }
    10911094      break;
    … …  
    10951098        print_double (NAN);
    10961099      } else {
    1097         print_double_exp (secfilt[seq].dMkron);
     1100        print_double_exp (secfilt[seq].dMkronChp);
    10981101      }
    10991102      break;

  • trunk/Ohana/src/opihi/dvo/imdata.c

    r39457 r40291  
    183183        for (i = 0; i < catalog.Nmeasure; i++) {
    184184          if ((catalog.measure[i].t < start) || (catalog.measure[i].t > stop)) continue;
    185           vec[0].elements.Flt[N] = catalog.measure[i].Mcal;
     185          vec[0].elements.Flt[N] = catalog.measure[i].McalPSF;
    186186          N++;
    187187          CHECK_REALLOCATE (vec[0].elements.Flt, opihi_flt, NPTS, N, 1000);

  • trunk/Ohana/src/opihi/dvo/imlist.c

    r39310 r40291  
    141141    if (VERBOSE) {
    142142      gprint (GP_LOG, "%3lld %s %8lld %8.4f %8.4f %f %5d %2d %4.2f %5.3f %5.3f",
    143                          (long long) i, image[i].name, (long long) image[i].imageID, r, d, t, image[i].nstar, image[i].photcode, image[i].secz, image[i].Mcal, image[i].dMcal);
     143                         (long long) i, image[i].name, (long long) image[i].imageID, r, d, t, image[i].nstar, image[i].photcode, image[i].secz, image[i].McalPSF, image[i].dMcal);
    144144
    145145      if (showUR) {

  • trunk/Ohana/src/opihi/dvo/imphot.c

    r39233 r40291  
    6464      for (x = 0; x < 100; x+=1.0, p++) {
    6565        // *p = applyMcal (&image[subset[0]], (fx*x), (fy*y));
    66         *p = image[subset[0]].Mcal;
     66        *p = image[subset[0]].McalPSF;
    6767      }
    6868    }
    … …  
    7171  for (j = 0; j < Nsubset; j++) {
    7272    i = subset[j];
    73     gprint (GP_ERR, "%s: %f\n", image[i].name, image[i].Mcal);
     73    gprint (GP_ERR, "%s: %f\n", image[i].name, image[i].McalPSF);
    7474
    7575// XXX old code when we had the option of a 2D zero point model
    … …  
    7777    switch (image[i].order) {
    7878    case 0:
    79       gprint (GP_ERR, "%s: %d - %f\n", image[i].name, image[i].order, image[i].Mcal);
     79      gprint (GP_ERR, "%s: %d - %f\n", image[i].name, image[i].order, image[i].McalPSF);
    8080      break;
    8181    case 1:
    82       gprint (GP_ERR, "%s: %d - %f, %d %d\n", image[i].name, image[i].order, image[i].Mcal, image[i].Mx, image[i].My);
     82      gprint (GP_ERR, "%s: %d - %f, %d %d\n", image[i].name, image[i].order, image[i].McalPSF, image[i].Mx, image[i].My);
    8383      break;
    8484    case 2:
    85       gprint (GP_ERR, "%s: %d - %f, %d %d, %d %d %d\n", image[i].name, image[i].order, image[i].Mcal, image[i].Mx, image[i].My, image[i].Mxx, image[i].Mxy, image[i].Myy);
     85      gprint (GP_ERR, "%s: %d - %f, %d %d, %d %d %d\n", image[i].name, image[i].order, image[i].McalPSF, image[i].Mx, image[i].My, image[i].Mxx, image[i].Mxy, image[i].Myy);
    8686      break;
    8787    case 3:
    88       gprint (GP_ERR, "%s: %d - %f, %d %d, %d %d %d, %d %d %d %d\n", image[i].name, image[i].order, image[i].Mcal, image[i].Mx, image[i].My,
     88      gprint (GP_ERR, "%s: %d - %f, %d %d, %d %d %d, %d %d %d %d\n", image[i].name, image[i].order, image[i].McalPSF, image[i].Mx, image[i].My,
    8989               image[i].Mxx, image[i].Mxy, image[i].Myy, image[i].Mxxx, image[i].Mxxy, image[i].Mxyy, image[i].Myyy);
    9090      break;
    9191    case 4:
    92       gprint (GP_ERR, "%s: %d - %f, %d %d, %d %d %d, %d %d %d %d, %d %d %d %d %d\n", image[i].name, image[i].order, image[i].Mcal, image[i].Mx, image[i].My,
     92      gprint (GP_ERR, "%s: %d - %f, %d %d, %d %d %d, %d %d %d %d, %d %d %d %d %d\n", image[i].name, image[i].order, image[i].McalPSF, image[i].Mx, image[i].My,
    9393               image[i].Mxx, image[i].Mxy, image[i].Myy, image[i].Mxxx, image[i].Mxxy, image[i].Mxyy, image[i].Myyy,
    9494               image[i].Mxxxx, image[i].Mxxxy, image[i].Mxxyy, image[i].Mxyyy, image[i].Myyyy);

  • trunk/Ohana/src/opihi/dvo/imstats.c

    r40165 r40291  
    4141    Xvec.elements.Flt[i] = image[i].secz;
    4242    if (Mcal)
    43       Yvec.elements.Flt[i] = image[i].Mcal;
     43      Yvec.elements.Flt[i] = image[i].McalPSF;
    4444    else
    4545      Yvec.elements.Flt[i] = image[i].dMcal;
    … …  
    4747    gprint (GP_ERR, "%d %8.4f %8.4f %10d %6d  %5.3f %6.3f %6.3f\n",
    4848             i, r, d, image[i].tzero, image[i].nstar, Xvec.elements.Flt[i],
    49              image[i].Mcal, image[i].dMcal);
     49             image[i].McalPSF, image[i].dMcal);
    5050  }
    5151  if (AutoLimits) SetLimits (&Xvec, &Yvec, &graphmode);

  • trunk/Ohana/src/opihi/dvo/objectcoverage.c

    r39457 r40291  
    203203        if (catalog.secfilt[j*Nsecfilt+Nsec].Ncode < 2) { continue; }
    204204
    205         invalid = ((catalog.secfilt[j*Nsecfilt + Nsec].M < 1.0) || (isnan(catalog.secfilt[j*Nsecfilt + Nsec].M)));
     205        invalid = ((catalog.secfilt[j*Nsecfilt + Nsec].MpsfChp < 1.0) || (isnan(catalog.secfilt[j*Nsecfilt + Nsec].MpsfChp)));
    206206        if (invalid) continue;
    207207       

  • trunk/Ohana/src/opihi/dvo/paverage.c

    r40165 r40291  
    125125      while (average[i].R > Rmax) average[i].R -= 360.0;
    126126
    127       mag = secfilt[i*Nsecfilt+Nsec].M;
     127      mag = secfilt[i*Nsecfilt+Nsec].MpsfChp;
    128128      Zvec[Npts] = MIN (1.0, MAX (0.01, (mag - Mz) / Mr));
    129129      if (LimExclude && (Zvec[Npts] > 0.99)) continue;

  • trunk/Ohana/src/opihi/dvo/remote.c

    r39283 r40291  
    2727    gprint (GP_ERR, "  -skip-result : do not try to read from the result file\n");
    2828    gprint (GP_ERR, "OR:    remote -reload (uniquer)\n");
     29    gprint (GP_ERR, "       (reloads the remote host results into vectors as if a parallel command were run)\n");
    2930    gprint (GP_ERR, "OR:    remote -get-results (uniquer)\n");
     31    gprint (GP_ERR, "       (generates the list of remote result filenames and status variables)\n");
     32    gprint (GP_ERR, "       (RESULT_FILE:i is the filenme, RESULT_STATUS:i is the dvo_client exit status)\n");
    3033    return FALSE;
    3134  }

  • trunk/Ohana/src/opihi/dvo/skycoverage.c

    r39233 r40291  
    330330              break;
    331331            case MIN_MCAL:
    332               V[ys*Nx + xs] = MIN(V[ys*Nx + xs], image[i].Mcal);
     332              V[ys*Nx + xs] = MIN(V[ys*Nx + xs], image[i].McalPSF);
    333333              break;
    334334            case MAX_MCAL:
    335               V[ys*Nx + xs] = MAX(V[ys*Nx + xs], image[i].Mcal);
     335              V[ys*Nx + xs] = MAX(V[ys*Nx + xs], image[i].McalPSF);
    336336              break;
    337337            case MIN_TIME: {

  • trunk/Ohana/src/opihi/include/astro.h

    r39610 r40291  
    4747double VectorFractionInterpolate (double *values, float fraction, int Npts);
    4848
    49 int PlxSetMeanEpoch (double *R, double *D, double *T, double *Rmean, double *Dmean, double *Tmean, int *mask, int Ntotal);
    50 int PlxSetEpochPosition (PlxFitData *fitdata, double *R, double *D, double *dR, double *dD, double *T, int *mask, int Ntotal, Coords *coords, double Tmean);
    51 int PlxOutlierClip (PlxFitData *fitdata, int *mask, int Noutlier, float dPsigMax, Vector *dPvec, int VERBOSE);
     49int PlxSetMeanEpoch (double *R, double *D, double *T, double *Rmean, double *Dmean, double *Tmean, opihi_int *mask, int Ntotal);
     50int PlxSetEpochPosition (PlxFitData *fitdata, double *R, double *D, double *dR, double *dD, double *T, opihi_int *mask, int Ntotal, Coords *coords, double Tmean);
     51int PlxOutlierClip (PlxFitData *fitdata, opihi_int *mask, int Noutlier, float dPsigMax, Vector *dPvec, int VERBOSE);
    5252
    5353int PlxFitDataAlloc (PlxFitData *data, int N);

  • trunk/Ohana/src/opihi/include/data.h

    r37807 r40291  
    162162/* starfuncs.c */
    163163double get_aperture_stats (Matrix *matrix, int X, int Y, int Npix, int Nborder, double max, int VERBOSE);
     164double get_box_stats (Matrix *matrix, int X, int Y, int dX, int dY, int Nborder, double max, int VERBOSE);
     165
    164166int set_rough_radii (double Ra, double Ri, double Ro);
    165167int get_rough_star (float *data, int Nx, int Ny, int x, int y, opihi_flt *xc, opihi_flt *yc, opihi_flt *sx, opihi_flt *sy, opihi_flt *sxy, opihi_flt *zs, opihi_flt *zp, opihi_flt *sk);

  • trunk/Ohana/src/opihi/lib.data/graphtools.c

    r38153 r40291  
    1010  if (xvec != NULL) {
    1111    if (xvec->type == OPIHI_FLT) {
    12       maxX = DBL_MIN;
     12      maxX = -DBL_MAX;
    1313      minX = DBL_MAX;
    1414      for (i = 0; i < xvec[0].Nelements; i++) {
    … …  
    3333  if (yvec != NULL) {
    3434    if (yvec->type == OPIHI_FLT) {
    35       maxY = DBL_MIN;
     35      maxY = -DBL_MAX;
    3636      minY = DBL_MAX;
    3737      for (i = 0; i < yvec[0].Nelements; i++) {

  • trunk/Ohana/src/opihi/lib.data/starfuncs.c

    r36679 r40291  
    101101}
    102102
     103double get_box_stats (Matrix *matrix, int X, int Y, int dX, int dY, int Nborder, double max, int VERBOSE) {
     104
     105  double *ring;
     106  double x, y, x2, y2, xy, I, sky, FWHMx, FWHMy, value, mag, Sxy;
     107  int i, j, n, Nring, Nmax;
     108  double Npts, gain, dsky2, dmag, peak, offset;
     109  char *string;
     110 
     111  string = get_variable ("GAIN");
     112  if (string == (char *) NULL) {
     113    gprint (GP_ERR, "assuming a value of 1.0\n");
     114    gain = 1.0;
     115  } else {
     116    gain = atof (string);
     117  }
     118  Nborder = MAX (1, Nborder);
     119  Nborder = MIN (1000, Nborder);
     120 
     121  int dX2 = (int)(0.5*dX);
     122  int dY2 = (int)(0.5*dY);
     123  dX = 2 * dX2 + 1;
     124  dY = 2 * dY2 + 1;
     125
     126  Nring = 2*Nborder*(dX + 2*Nborder) + 2*Nborder*(dY + 2*Nborder);
     127  ALLOCATE (ring, double, Nring);
     128  bzero (ring, sizeof(double)*Nring);
     129
     130  // get the pixels in the border regions:
     131  // XXX gfits_get_matrix_value returns 0 for out-of-bounds pixels, but should return NAN
     132  // and they should be skipped
     133  n = 0; 
     134  for (j = 0; j < Nborder; j++) {
     135    for (i = X - dX2 - Nborder; i < X + dX2 + Nborder + 1; i++) {
     136      value = gfits_get_matrix_value (matrix, i, (int)(Y - dY2 - j));
     137      if (isfinite(value)) { ring[n] = value; n++; }
     138      value = gfits_get_matrix_value (matrix, i, (int)(Y + dY2 + j));
     139      if (isfinite(value)) { ring[n] = value; n++; }
     140    }
     141    for (i = Y - dY2; i < Y + dY2 + 1; i++) {
     142      value = gfits_get_matrix_value (matrix, (int)(X - dX2 - j), i);
     143      if (isfinite(value)) { ring[n] = value; n++; }
     144      value = gfits_get_matrix_value (matrix, (int)(X + dX2 + j), i);
     145      if (isfinite(value)) { ring[n] = value; n++; }
     146    }
     147  }
     148  Nring = n;
     149  dsort (ring, Nring);
     150  for (Npts = sky = dsky2 = 0, i = 0.25*Nring; i < 0.75*Nring; i++, Npts += 1.0) {
     151    sky += ring[i];
     152    dsky2 += ring[i]*ring[i];
     153  }
     154  sky = sky / Npts;
     155  dsky2 = dsky2 / Npts - sky*sky;
     156  free (ring);
     157
     158  float dx, dy;
     159
     160  peak = 0;
     161  Npts = Nmax = 0;
     162  x = y = x2 = y2 = xy = I = 0;
     163  for (i = X - dX2; i < X + dX2 + 1; i++) {
     164    for (j = Y - dY2; j < Y + dY2 + 1; j++) {
     165      value = gfits_get_matrix_value (matrix, i, j);
     166      if (!isfinite(value)) continue;
     167      offset = value - sky;
     168      dx = i - X;
     169      dy = j - Y;
     170      x  += dx*offset;
     171      y  += dy*offset;
     172      x2 += dx*dx*offset;
     173      y2 += dy*dy*offset;
     174      xy += dx*dy*offset;
     175      I  += offset;
     176      Npts ++;
     177      if (value > max) {
     178        Nmax ++;
     179      }
     180      if (value > peak) peak = value;
     181    }
     182  }
     183
     184  x = x / I;
     185  y = y / I;
     186  FWHMx = 2.355*sqrt (fabs(x2 / I - x*x));
     187  FWHMy = 2.355*sqrt (fabs(y2 / I - y*y));
     188  Sxy   = xy / I - x*y;
     189  mag = -2.5*log10(I);
     190
     191  // flux_error = sqrt( I + Npts*dsky2 )
     192  // dmag = 1.086 * flux_error / flux
     193  dmag = 1.086 * sqrt (fabs(I + Npts*dsky2)) / (gain * I);
     194  x = x + X;
     195  y = y + Y;
     196 
     197  set_variable ("Xg", x);
     198  set_variable ("Yg", y);
     199  set_variable ("SXg", FWHMx);
     200  set_variable ("SYg", FWHMy);
     201  set_variable ("SXYg", Sxy);
     202  set_variable ("Sg", sky);
     203  set_variable ("dSg", sqrt (fabs (dsky2)));
     204  set_variable ("Zg", mag);
     205  set_variable ("dZg", dmag);
     206  set_variable ("Zcg", I);
     207  set_variable ("Zpk", peak);
     208  set_int_variable ("Nsat", Nmax);
     209  set_int_variable ("Npts", Npts);
     210 
     211  if (VERBOSE) gprint (GP_LOG, "%f %f %f %f %f %f %f %f\n", x, y, FWHMx, FWHMy, sky, I, mag, dmag);
     212
     213  return (mag);
     214
     215}
     216
    103217static double Raper  =  5;
    104218static double Rinner = 10;

  • trunk/Ohana/src/opihi/lib.shell/VectorIO.c

    r39457 r40291  
    4141    for (j = 0; j < Nvec; j++) {
    4242      // if the format is not defined, just use the native byte-widths
    43       tformat[2*j + 0] = (vec[j][0].type == OPIHI_FLT) ? 'D' : 'J';
     43      tformat[2*j + 0] = (vec[j][0].type == OPIHI_FLT) ? 'D' : 'K'; // this depends on opihi_int == int64_t for Int
    4444      tformat[2*j + 1] = 0;
    4545    }
    … …  
    6060  for (j = 0; j < Nvec; j++) {
    6161    if (vec[j][0].type == OPIHI_FLT) {
    62       gfits_set_bintable_column_reformat (theader, ftable, vec[j][0].name, "double", vec[j][0].elements.Flt, vec[j][0].Nelements, nativeOrder);
     62      gfits_set_bintable_column_reformat (theader, ftable, vec[j][0].name, "double",  vec[j][0].elements.Flt, vec[j][0].Nelements, nativeOrder);
    6363    } else {
    64       gfits_set_bintable_column_reformat (theader, ftable, vec[j][0].name, "int", vec[j][0].elements.Int, vec[j][0].Nelements, nativeOrder);
     64//    gfits_set_bintable_column_reformat (theader, ftable, vec[j][0].name, "int",     vec[j][0].elements.Int, vec[j][0].Nelements, nativeOrder);
     65      gfits_set_bintable_column_reformat (theader, ftable, vec[j][0].name, "int64_t", vec[j][0].elements.Int, vec[j][0].Nelements, nativeOrder);
    6566    }
    6667  }
    … …  
    328329  ASSIGN_DATA(short,   short,   Int);
    329330  ASSIGN_DATA(int,     int,     Int);
    330   ASSIGN_DATA(int64_t, int64_t, Flt); // int64_t has a problem: Int is too small, Flt is wrong precision
     331  ASSIGN_DATA(int64_t, int64_t, Int); // XXX this works if opihi_int is assigned to int64_t
     332//ASSIGN_DATA(int64_t, int64_t, Flt); // int64_t has a problem: Int is too small, Flt is wrong precision
    331333  ASSIGN_DATA(float,   float,   Flt);
    332334  ASSIGN_DATA(double,  double,  Flt);
    … …  
    353355  ASSIGN_DATA_TRANSPOSE(short,   short,   Int);
    354356  ASSIGN_DATA_TRANSPOSE(int,     int,     Int);
    355   ASSIGN_DATA_TRANSPOSE(int64_t, int64_t, Flt);
     357  ASSIGN_DATA_TRANSPOSE(int64_t, int64_t, Int);
     358//ASSIGN_DATA_TRANSPOSE(int64_t, int64_t, Flt); // see above comment
    356359  ASSIGN_DATA_TRANSPOSE(float,   float,   Flt);
    357360  ASSIGN_DATA_TRANSPOSE(double,  double,  Flt);

  • trunk/Ohana/src/opihi/lib.shell/convert_to_RPN.c

    r39558 r40291  
    120120        Nop_stack ++;
    121121        break;
    122       case ST_UNARY:
    123122      case ST_BINARY:
    124123      case ST_TRINARY:
    … …  
    139138        Nop_stack ++;
    140139        break;
     140      case ST_UNARY:
    141141      case ST_LEFT: 
    142142        /* push operator on OP stack */

  • trunk/Ohana/src/opihi/lib.shell/dvomath.c

    r39457 r40291  
    2929  unsigned int Ncstack;
    3030  cstack = isolate_elements (argc, argv, &Ncstack);
     31
     32  // for (i = 0; i < Ncstack; i++) {
     33  //   fprintf (stderr, "%d : %s\n", i, cstack[i]);
     34  // }
    3135
    3236  /* generate RPN stack from cstack arguments */
    … …  
    8286      } else {
    8387        if (stack[0].type == ST_SCALAR_INT) {
    84           sprintf (outname, "%d", stack[0].IntValue);
     88          sprintf (outname, OPIHI_INT_FMT, stack[0].IntValue);
    8589        } else {
    8690          sprintf (outname, "%.12g", stack[0].FltValue);

  • trunk/Ohana/src/opihi/lib.shell/evaluate_stack.c

    r40014 r40291  
    8989      }
    9090      if (tmp_stack.type == ST_SCALAR_INT) {
    91         gprint (GP_ERR, "---> %d ", tmp_stack.IntValue);
     91        gprint (GP_ERR, "---> "OPIHI_INT_FMT" ", tmp_stack.IntValue);
    9292      }
    9393      if (tmp_stack.type == ST_SCALAR_FLT) {
    … …  
    110110
    111111        if (i < 3) {  /* need two variables to operate on */
    112           snprintf (line, 512, "syntax error: trinary operator without three operands: %s\n", stack[i].name);
     112          snprintf (line, 512, "syntax error: trinary operator without three operands: %s\n(Note that the : in a trinary operation must be protected by spaces", stack[i].name);
    113113          push_error (line);
    114114          clear_stack (&tmp_stack);
    … …  
    124124
    125125        /* there are no valid unary string operators */
    126         snprintf (line, 512, "invalid operands for trinary operator %s (mismatch types?)", stack[i].name);
     126        snprintf (line, 512, "invalid operands for trinary operator %s (mismatch types?)\n(Note that the : in a trinary operation must be protected by spaces)", stack[i].name);
    127127        push_error (line);
    128128        clear_stack (&tmp_stack);
    … …  
    131131      got_three_op:
    132132        if (!status) {
    133           snprintf (line, 512, "syntax error: invalid operand for trinary operation: %s or %s or %s\n", stack[i-1].name, stack[i-2].name, stack[i-3].name);
     133          snprintf (line, 512, "syntax error: invalid operand for trinary operation: %s or %s or %s\n(Note that the : in a trinary operation must be protected by spaces)", stack[i-1].name, stack[i-2].name, stack[i-3].name);
    134134          push_error (line);
    135135          clear_stack (&tmp_stack);
    … …  
    188188          snprintf (line, 512, "syntax error: invalid operand for binary operation: %s or %s\n", stack[i-1].name, stack[i-2].name);
    189189          push_error (line);
     190          if (strchr(stack[i-1].name, ':') || strchr(stack[i-2].name, ':')) {
     191            snprintf (line, 512, "syntax error: invalid operand for binary operation: %s or %s\n(Note that the : in a trinary operation must be protected by spaces)\n", stack[i-1].name, stack[i-2].name);
     192            push_error (line);
     193          }
    190194          clear_stack (&tmp_stack);
    191195          return (FALSE);

  • trunk/Ohana/src/opihi/lib.shell/parse.c

    r33662 r40291  
    245245        vec[0].elements.Flt[Nx] = atof (val);
    246246      } else {
    247         vec[0].elements.Int[Nx] = atol (val);
     247        vec[0].elements.Int[Nx] = atoll (val);
    248248      }
    249249    }
Note: See TracChangeset for help on using the changeset viewer.