IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 39612 for trunk


Ignore:
Timestamp:
Jun 25, 2016, 11:45:39 AM (10 years ago)
Author:
eugene
Message:

updated FitPM, FitPosPMfixed, FitPMandPar to use the modified concept for IRLS: the weight factor multiplies the error

Location:
trunk/Ohana/src/relastro
Files:
1 added
8 edited

Legend:

Unmodified
Added
Removed

  • trunk/Ohana/src/relastro/Makefile

    r39474 r39612  
    1 default: relastro relastro_client fitobj fitobj2 fitpm
     1default: relastro relastro_client fitobj fitobj2 fitobj3 fitpm
    22help:
    33        @echo "make options: relastro (default)"
    … …  
    238238$(BIN)/fitobj2.$(ARCH): $(FITOBJ2)
    239239
     240fitobj3: $(BIN)/fitobj3.$(ARCH)
     241
     242FITOBJ3 = \
     243$(SRC)/fitobj3.$(ARCH).o              \
     244$(SRC)/ParFactor.$(ARCH).o           \
     245$(SRC)/FitAstromOps.$(ARCH).o        \
     246$(SRC)/FitPosPMfixed.$(ARCH).o               \
     247$(SRC)/FitPM.$(ARCH).o               \
     248$(SRC)/FitPMandPar.$(ARCH).o               \
     249$(SRC)/FitPosStack.$(ARCH).o       \
     250$(SRC)/UpdateObjects.$(ARCH).o               \
     251$(SRC)/BootstrapOps.$(ARCH).o        \
     252$(SRC)/mkpolyterm.$(ARCH).o            \
     253$(SRC)/GetAstromError.$(ARCH).o            \
     254$(SRC)/MeasFilterTest.$(ARCH).o            \
     255$(SRC)/extra.$(ARCH).o            \
     256$(SRC)/plotstuff.$(ARCH).o            \
     257$(SRC)/fitpoly.$(ARCH).o
     258
     259$(FITOBJ3): $(INC)/relastro.h
     260$(BIN)/fitobj3.$(ARCH): $(FITOBJ3)
     261
    240262fitstk: $(BIN)/fitstk.$(ARCH)
    241263

  • trunk/Ohana/src/relastro/include/relastro.h

    r39610 r39612  
    55# include <assert.h>
    66# include <pthread.h>
     7
     8# define OLD_METHOD 0
    79
    810// choose off_t or int depending on full-scale relphot analysis resources

  • trunk/Ohana/src/relastro/src/FitPM.c

    r39371 r39612  
    11# include "relastro.h"
     2
    23int FitPM_MinChisq (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints);
    34int FitPM_SetChisq (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints);
    … …  
    1617  // Convert the measurement errors into initial weights.
    1718  for (i = 0; i < Npoints; i++) {
    18     points[i].Wx = 1.0 / SQ(points[i].dX);
    19     points[i].Wy = 1.0 / SQ(points[i].dY);
     19    points[i].Wx = 1.0;
     20    points[i].Wy = 1.0;
     21    points[i].Qx = 1.0 / SQ(points[i].dX);
     22    points[i].Qy = 1.0 / SQ(points[i].dY);
     23    points[i].qx = points[i].Wx * points[i].Qx; // Wx, Wy start out at 1.0
     24    points[i].qy = points[i].Wy * points[i].Qy; // Wx, Wy start out at 1.0
    2025  }
    2126
    … …  
    3338  // Convert the measurement errors into initial weights.
    3439  for (i = 0; i < Npoints; i++) {
    35     points[i].Wx = 1.0 / SQ(points[i].dX);
    36     points[i].Wy = 1.0 / SQ(points[i].dY);
     40    points[i].Wx = 1.0;
     41    points[i].Wy = 1.0;
     42    points[i].Qx = 1.0 / SQ(points[i].dX);
     43    points[i].Qy = 1.0 / SQ(points[i].dY);
     44    points[i].qx = points[i].Wx * points[i].Qx; // Wx, Wy start out at 1.0
     45    points[i].qy = points[i].Wy * points[i].Qy; // Wx, Wy start out at 1.0
    3746  }
    3847 
    … …  
    7685      points[i].Wx = weight_cauchy(points[i].rx / points[i].dX);
    7786      points[i].Wy = weight_cauchy(points[i].ry / points[i].dY);
     87      points[i].qx = points[i].Wx * points[i].Qx;
     88      points[i].qy = points[i].Wy * points[i].Qy;
    7889    }   
    7990
    … …  
    201212    Nfit ++;
    202213
    203     wx = points[i].Wx;
    204     wy = points[i].Wy;
     214    wx = points[i].qx;
     215    wy = points[i].qy;
    205216
    206217    Wx += wx;
    … …  
    298309  // the reduced chisq is divided by (Ndof = 2*Nfit - Nterms)
    299310  fit[0].chisq = chisq / (2.0*fit[0].Nfit - data->Nterms);
    300 
    301311  return TRUE;
    302312}

  • trunk/Ohana/src/relastro/src/FitPMandPar.c

    r39610 r39612  
    11# include "relastro.h"
    2 # define OLD_METHOD 0
    32
    43int FitPMandPar_MinChisq (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints);
    … …  
    1817  // Convert the measurement errors into initial weights.
    1918  for (i = 0; i < Npoints; i++) {
    20 # if (OLD_METHOD)
    21     points[i].Wx = 1.0 / SQ(points[i].dX);
    22     points[i].Wy = 1.0 / SQ(points[i].dY);
    23 # else
    2419    points[i].Wx = 1.0;
    2520    points[i].Wy = 1.0;
    … …  
    2823    points[i].qx = points[i].Wx * points[i].Qx; // Wx, Wy start out at 1.0
    2924    points[i].qy = points[i].Wy * points[i].Qy; // Wx, Wy start out at 1.0
    30 # endif
    3125  }
    3226
    3327  if (!FitPMandPar_MinChisq (fit, data, points, Npoints)) return FALSE;
    34 
    3528  if (!FitPMandPar_SetChisq (fit, data, points, Npoints)) return FALSE;
    3629  return TRUE;
    … …  
    4538  // Convert the measurement errors into initial weights.
    4639  for (i = 0; i < Npoints; i++) {
    47 # if (OLD_METHOD)
    48     points[i].Wx = 1.0 / SQ(points[i].dX);
    49     points[i].Wy = 1.0 / SQ(points[i].dY);
    50     points[i].qx = points[i].Wx; // Wx, Wy start out at 1.0
    51     points[i].qy = points[i].Wy; // Wx, Wy start out at 1.0
    52 # else
     40    points[i].Wx = 1.0;
     41    points[i].Wy = 1.0;
    5342    points[i].Qx = 1.0 / SQ(points[i].dX);
    5443    points[i].Qy = 1.0 / SQ(points[i].dY);
    5544    points[i].qx = points[i].Wx * points[i].Qx; // Wx, Wy start out at 1.0
    5645    points[i].qy = points[i].Wy * points[i].Qy; // Wx, Wy start out at 1.0
    57 # endif
    5846  }
    5947 
    … …  
    231219    Nfit ++;
    232220
    233 # if (OLD_METHOD)   
    234     wx = points[i].Wx;
    235     wy = points[i].Wy;
    236 # else
    237221    wx = points[i].qx;
    238222    wy = points[i].qy;
    239 # endif
    240223
    241224    Wx += wx;

  • trunk/Ohana/src/relastro/src/FitPosPMfixed.c

    r39375 r39612  
    11# include "relastro.h"
     2
    23int FitPosPMfixed_MinChisq (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints);
    34int FitPosPMfixed_SetChisq (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints);
    … …  
    4546  // Convert the measurement errors into initial weights.
    4647  for (i = 0; i < Npoints; i++) {
    47     points[i].Wx = 1.0 / SQ(points[i].dX);
    48     points[i].Wy = 1.0 / SQ(points[i].dY);
     48    points[i].Wx = 1.0;
     49    points[i].Wy = 1.0;
     50    points[i].Qx = 1.0 / SQ(points[i].dX);
     51    points[i].Qy = 1.0 / SQ(points[i].dY);
     52    points[i].qx = points[i].Wx * points[i].Qx; // Wx, Wy start out at 1.0
     53    points[i].qy = points[i].Wy * points[i].Qy; // Wx, Wy start out at 1.0
    4954  }
    5055
    … …  
    6974  // Convert the measurement errors into initial weights.
    7075  for (i = 0; i < Npoints; i++) {
    71     points[i].Wx = 1.0 / SQ(points[i].dX);
    72     points[i].Wy = 1.0 / SQ(points[i].dY);
     76    points[i].Wx = 1.0;
     77    points[i].Wy = 1.0;
     78    points[i].Qx = 1.0 / SQ(points[i].dX);
     79    points[i].Qy = 1.0 / SQ(points[i].dY);
     80    points[i].qx = points[i].Wx * points[i].Qx; // Wx, Wy start out at 1.0
     81    points[i].qy = points[i].Wy * points[i].Qy; // Wx, Wy start out at 1.0
    7382  }
    7483 
    … …  
    112121      points[i].Wx = weight_cauchy(points[i].rx / points[i].dX);
    113122      points[i].Wy = weight_cauchy(points[i].ry / points[i].dY);
     123      points[i].qx = points[i].Wx * points[i].Qx;
     124      points[i].qy = points[i].Wy * points[i].Qy;
    114125    }   
    115126
    … …  
    231242    Nfit ++;
    232243
    233     wx = points[i].Wx;
    234     wy = points[i].Wy;
     244    wx = points[i].qx;
     245    wy = points[i].qy;
    235246
    236247    Wx += wx;

  • trunk/Ohana/src/relastro/src/fitobj.c

    r39610 r39612  
    6969  if (argc != 2) {
    7070    fprintf (stderr, "USAGE: %s (mode)\n", argv[0]);
    71     fprintf (stderr, "  mode: pos, pos-irls, pos-boot, pm, pm-irls, pm-boot, plx, plx-irls, plx-boot\n");
     71    fprintf (stderr, "  mode: pos, pm, plx\n");
     72    fprintf (stderr, "  test astrometry fitting using UpdateObjects_Chip code\n");
    7273    fprintf (stderr, " options: -Nstars [3000], -Npoints [100], -Noutliers [10], -Nbootstrap [100]\n");
    7374    exit (2);
    … …  
    9899  // load reference table here
    99100  if (!LoadPhotcodesText ("dvo.photcodes")) {
    100     fprintf (stderr, "failed to load photcodes\n");
     101    fprintf (stderr, "failed to load photcode file dvo.photcodes (copy from ipp/ippconfig)\n");
    101102    exit (2);
    102103  }
    … …  
    121122  for (i = 0; i < Nstars; i++) {
    122123
     124    // generate random stars uniformly distributed on the sphere (ra = 0,36)
    123125    double Ro  = 360.0*drand48();
    124126    double Phi = 2.0*(drand48() - 0.5);

  • trunk/Ohana/src/relastro/src/fitobj2.c

    r39457 r39612  
    4848  int N_STARS     = 3000;
    4949  int N_POINTS    =  100;
     50  int N_MIN_POINTS =  10;
    5051  int N_BOOTSTRAP =  100;
    5152  float F_OUTLIERS  =  0.1;
    … …  
    6263    remove_argument (N, &argc, argv);
    6364  }
     65  if ((N = get_argument (argc, argv, "-Nminpoints"))) {
     66    remove_argument (N, &argc, argv);
     67    N_MIN_POINTS = atoi (argv[N]);
     68    remove_argument (N, &argc, argv);
     69  }
    6470  if ((N = get_argument (argc, argv, "-foutliers"))) {
    6571    remove_argument (N, &argc, argv);
    … …  
    7682  if (argc != 2) {
    7783    fprintf (stderr, "USAGE: %s (mode)\n", argv[0]);
    78     fprintf (stderr, "  mode: pos, pos-irls, pos-boot, pm, pm-irls, pm-boot, plx, plx-irls, plx-boot\n");
    79     fprintf (stderr, " options: -Nstars [3000], -Npoints [100], -Noutliers [10], -Nbootstrap [100]\n");
     84    fprintf (stderr, "  mode: pos, pm, plx\n");
     85    fprintf (stderr, "  test astrometry fitting using UpdateObjects_Chip code (random number of measurements)\n");
     86    fprintf (stderr, " options: -Nstars [3000], -Npoints [100], -foutliers [0.1], -Nbootstrap [100]\n");
    8087    exit (2);
    8188  }
    … …  
    94101    fprintf (stderr, "USAGE: %s (mode)\n", argv[0]);
    95102    fprintf (stderr, "  mode: pos, pm, plx\n");
    96     fprintf (stderr, " options: -Nstars [3000], -Npoints [100], -Noutliers [10], -Nbootstrap [100]\n");
     103    fprintf (stderr, " options: -Nstars [3000], -Npoints [100], -foutliers [0.1], -Nbootstrap [100]\n");
    97104    exit (2);
    98105  }
    … …  
    105112  // load reference table here
    106113  if (!LoadPhotcodesText ("dvo.photcodes")) {
    107     fprintf (stderr, "failed to load photcodes\n");
     114    fprintf (stderr, "failed to load photcode file dvo.photcodes (copy from ipp/ippconfig)\n");
    108115    exit (2);
    109116  }
    … …  
    128135  for (i = 0; i < Nstars; i++) {
    129136
    130     double Ro  = 360.0*(drand48_cnt() - 0.5);
     137    // generate random stars uniformly distributed on the sphere (ra = 0,36)
     138    double Ro  = 360.0*drand48_cnt();
    131139    double Phi = 2.0*(drand48_cnt() - 0.5);
    132140    double Do  = DEG_RAD * asin(Phi);
    … …  
    138146    if (FIT_MODE != FIT_PM_AND_PAR) plx = 0.0;
    139147
    140     int Npoints = N_POINTS*drand48_cnt(); // minimum of 0 points
     148    int Npoints = (N_POINTS - N_MIN_POINTS)*drand48_cnt() + N_MIN_POINTS; // randomly choose N_MIN_POINTS to N_POINTS points
    141149    int Noutliers_max = Npoints * F_OUTLIERS;
    142150    int Noutliers = Noutliers_max*drand48_cnt(); // minimum of 0 points

  • trunk/Ohana/src/relastro/test/mana.sh

    r39610 r39612  
    4040end
    4141
     42macro testobj3
     43  if ($0 != 3)
     44    echo "USAGE: testobj3 (Nstars) (foutliers)"
     45    break
     46  end
     47
     48  $Nmin:pos = 99
     49  $Nmin:pm  = 99
     50  $Nmin:plx = 99
     51
     52  foreach mode pos pm plx
     53    exec ../bin/fitobj3.lin64 $mode -Nstars $1 -foutliers $2 -Nminpoints $Nmin:$mode > test.$mode.dat
     54    checkstats test.$mode.dat $mode
     55  end
     56end
     57
    4258macro testobj2
    4359  if ($0 != 3)
    … …  
    4662  end
    4763
    48   exec ../bin/fitobj2.lin64 plx -Nstars $1 -foutliers $2 > test.plx.dat
    49   checkstats test.plx.dat
     64  $Nmin:pos = 1
     65  $Nmin:pm  = 5
     66  $Nmin:plx = 7
     67
     68  foreach mode pos pm plx
     69    exec ../bin/fitobj2.lin64 $mode -Nstars $1 -foutliers $2 -Nminpoints $Nmin:$mode > test.$mode.dat
     70    checkstats test.$mode.dat $mode
     71  end
    5072end
    5173
    … …  
    5880  foreach mode pos pm plx
    5981    exec ../bin/fitobj.lin64 $mode > test.$mode.dat
    60     if ("$mode" == "plx")
    61       echo "is plx"
    62       reload test.$mode.dat 1
    63     else
    64       reload test.$mode.dat 0
    65     end
    66 
     82    checkstats.obj test.$mode.dat $mode
    6783    fprintf "%5.2f %5.2f %5.2f %5.2f %5.2f : %5.2f %5.2f %5.2f %5.2f %5.2f : %5.2f %5.2f %5.2f %5.2f %5.2f : %6.2f %5.2f" $value_dRf $value_dDf $value_duR $value_duD $value_dP $value_dRo $value_dDo $value_duRo $value_duDo $value_dPo $sigma_dRo $sigma_dDo $sigma_duRo $sigma_duDo $sigma_dPo $mean_Nfit $mean_chisq
    6884  end
    … …  
    108124  data $1
    109125  if ($2)
    110     # read Npts 1 Nout 2 Ro 3 Do 4 uRo 5 uDo 6 Po 7 Rx 9 Dx 10 uRx 11 uDx 12 Px 13 Tx 15 dRo 17 dDo 18 duRo 19 duDo 20 dPo 21 Nfit 22 chisq 26
     126    read Npts 1 Nout 2 Ro 3 Do 4 uRo 5 uDo 6 Po 7 Rx 9 Dx 10 uRx 11 uDx 12 Px 13 Tx 15 dRo 17 dDo 18 duRo 19 duDo 20 dPo 21 Nfit 22 chisq 26
    111127    read Ro 1 Do 2 uRo 3 uDo 4 Po 5 Rx 7 Dx 8 uRx 9 uDx 10 Px 11 Tx 13 dRo 15 dDo 16 duRo 17 duDo 18 dPo 19 Nfit 20 chisq 24
    112128  else
    … …  
    145161end
    146162
    147 macro teststats
    148   if ($0 != 2)
    149     echo "USAGE: teststats (file)"
     163macro checkstats.obj
     164  if ($0 != 3)
     165    echo "USAGE: checkstats.obj (file) (mode)"
    150166    break
    151167  end
    152168
    153169  data $1
    154   read Npts 1 Nout 2 Ro 3 Do 4 uRo 5 uDo 6 Po 7 Rx 9 Dx 10 uRx 11 uDx 12 Px 13 Tx 15 dRo 17 dDo 18 duRo 19 duDo 20 dPo 21 Nfit 22 chisq 26
    155   # read Ro 1 Do 2 uRo 3 uDo 4 Po 5 Rx 7 Dx 8 uRx 9 uDx 10 Px 11 Tx 13 dRo 15 dDo 16 duRo 17 duDo 18 dPo 19 Nfit 20 chisq 24
    156   # read Ro 1 Do 2 uRo 3 uDo 4 Po 5 Rx 7 Dx 8 uRx 9 uDx 10 Px 11 Tx 13 dRo 15 dDo 16 duRo 17 duDo 18 dPo 19 Nfit 20 chisq 22
     170  if ("$2" == "pos")
     171    read Ro 1 Do 2 uRo 3 uDo 4 Po 5 Rx 7 Dx 8 uRx 9 uDx 10 Px 11 Tx 13 dRo 15 dDo 16 duRo 17 duDo 18 dPo 19 Nfit 20 chisq 22
     172  end
     173  if ("$2" == "pm")
     174    read Ro 1 Do 2 uRo 3 uDo 4 Po 5 Rx 7 Dx 8 uRx 9 uDx 10 Px 11 Tx 13 dRo 15 dDo 16 duRo 17 duDo 18 dPo 19 Nfit 20 chisq 23
     175  end
     176  if ("$2" == "plx")
     177    read Ro 1 Do 2 uRo 3 uDo 4 Po 5 Rx 7 Dx 8 uRx 9 uDx 10 Px 11 Tx 13 dRo 15 dDo 16 duRo 17 duDo 18 dPo 19 Nfit 20 chisq 24
     178  end
    157179
    158180  foreach f R D uD uR P
    … …  
    185207  # echo chisq
    186208  vstat -q chisq; $mean_chisq = $MEAN
    187 
    188   fprintf "%5.2f %5.2f %5.2f %5.2f %5.2f : %5.2f %5.2f %5.2f %5.2f %5.2f : %5.2f %5.2f %5.2f %5.2f %5.2f : %6.2f %5.2f" $value_dRf $value_dDf $value_duR $value_duD $value_dP $value_dRo $value_dDo $value_duRo $value_duDo $value_dPo $sigma_dRo $sigma_dDo $sigma_duRo $sigma_duDo $sigma_dPo $mean_Nfit $mean_chisq
    189 
    190 end
    191 
    192 macro checkstats
    193   if ($0 != 2)
    194     echo "USAGE: checkstats (file)"
    195     break
    196   end
    197 
    198   teststats $1
    199  
    200   delete dRbin oRbin
    201   for i 0 100
    202     subset tmp = dRf if (Npts == $i)
    203     vstat tmp -q
    204     concat $SIGMA dRbin
    205     concat $MEAN oRbin
    206   end
    207 
    208   create n 0 dRbin[]
    209   set dRmod = 0.010 / sqrt (n)
    210   lim n -0.001 0.015; clear; box; plot n dRbin
    211   plot -c red -pt 7 n dRmod
    212 end
    213 
    214 macro load.stack
    215   if ($0 != 2)
    216     echo "USAGE: load.stack (file)"
     209end
     210
     211macro teststats
     212  if ($0 != 3)
     213    echo "USAGE: teststats (file) (mode)"
    217214    break
    218215  end
    219216
    220217  data $1
    221   read Npts 1 Nout 2 Nstk 3 Ro 4 Do 5 uRo 6 uDo 7 Po 8 Rx 10 Dx 11 uRx 12 uDx 3 Px 14 Tx 16 dRo 18 dDo 19 duRo 20 duDo 21 dPo 22 Nfit 23 Rstk 25 Dstk 26 dRstk 27 dDstk 28 chisq 32
    222   # read Ro 1 Do 2 uRo 3 uDo 4 Po 5 Rx 7 Dx 8 uRx 9 uDx 10 Px 11 Tx 13 dRo 15 dDo 16 duRo 17 duDo 18 dPo 19 Nfit 20 chisq 24
    223   # read Ro 1 Do 2 uRo 3 uDo 4 Po 5 Rx 7 Dx 8 uRx 9 uDx 10 Px 11 Tx 13 dRo 15 dDo 16 duRo 17 duDo 18 dPo 19 Nfit 20 chisq 22
    224 
    225   foreach f uD uR P
     218
     219  if ("$2" == "pos")
     220    read Npts 1 Nout 2 Ro 3 Do 4 uRo 5 uDo 6 Po 7 Rx 9 Dx 10 uRx 11 uDx 12 Px 13 Tx 15 dRo 17 dDo 18 duRo 19 duDo 20 dPo 21 Nfit 22 chisq 24
     221  end
     222  if ("$2" == "pm")
     223    read Npts 1 Nout 2 Ro 3 Do 4 uRo 5 uDo 6 Po 7 Rx 9 Dx 10 uRx 11 uDx 12 Px 13 Tx 15 dRo 17 dDo 18 duRo 19 duDo 20 dPo 21 Nfit 22 chisq 25
     224  end
     225  if ("$2" == "plx")
     226    read Npts 1 Nout 2 Ro 3 Do 4 uRo 5 uDo 6 Po 7 Rx 9 Dx 10 uRx 11 uDx 12 Px 13 Tx 15 dRo 17 dDo 18 duRo 19 duDo 20 dPo 21 Nfit 22 chisq 26
     227  end
     228
     229  foreach f R D uD uR P
    226230    set d$f = $f\o - $f\x
    227231  end
    228   set sdR = 3600*(Ro - Rstk)*dcos(Do)
    229   set sdD = 3600*(Do - Dstk)
    230232
    231233  set Rf = Rx - uRx*(Tx - $Tref)/3600.0/dcos(Dx)
    … …  
    242244  vstat -q dP ;  $value_dP  =  $SIGMA*1000
    243245
    244   vstat -q sdR;  $value_sdR =  $SIGMA*1000
    245   vstat -q sdD;  $value_sdD =  $SIGMA*1000
    246 
    247246  # echo "fit errors"
    248247  vstat -q dRo;  $value_dRo  = $MEAN*1000; $sigma_dRo  = $SIGMA*1000
    … …  
    259258
    260259  fprintf "%5.2f %5.2f %5.2f %5.2f %5.2f : %5.2f %5.2f %5.2f %5.2f %5.2f : %5.2f %5.2f %5.2f %5.2f %5.2f : %6.2f %5.2f" $value_dRf $value_dDf $value_duR $value_duD $value_dP $value_dRo $value_dDo $value_duRo $value_duDo $value_dPo $sigma_dRo $sigma_dDo $sigma_duRo $sigma_duDo $sigma_dPo $mean_Nfit $mean_chisq
     260
     261end
     262
     263macro checkstats
     264  if ($0 != 3)
     265    echo "USAGE: checkstats (file) (mode)"
     266    break
     267  end
     268
     269  teststats $1 $2
    261270 
    262   delete dRbin oRbin
     271  delete -q dRbin oRbin
     272  for i 0 100
     273    subset tmp = dRf if (Npts == $i)
     274    vstat tmp -q
     275    concat $SIGMA dRbin
     276    concat $MEAN oRbin
     277  end
     278
     279  create n 0 dRbin[]
     280  set dRmod = 0.010 / sqrt (n)
     281  lim n -0.001 0.015; clear; box; plot n dRbin
     282  plot -c red -pt 7 n dRmod
     283end
     284
     285macro load.stack
     286  if ($0 != 2)
     287    echo "USAGE: load.stack (file)"
     288    break
     289  end
     290
     291  data $1
     292  read Npts 1 Nout 2 Nstk 3 Ro 4 Do 5 uRo 6 uDo 7 Po 8 Rx 10 Dx 11 uRx 12 uDx 3 Px 14 Tx 16 dRo 18 dDo 19 duRo 20 duDo 21 dPo 22 Nfit 23 Rstk 25 Dstk 26 dRstk 27 dDstk 28 chisq 32
     293  # read Ro 1 Do 2 uRo 3 uDo 4 Po 5 Rx 7 Dx 8 uRx 9 uDx 10 Px 11 Tx 13 dRo 15 dDo 16 duRo 17 duDo 18 dPo 19 Nfit 20 chisq 24
     294  # read Ro 1 Do 2 uRo 3 uDo 4 Po 5 Rx 7 Dx 8 uRx 9 uDx 10 Px 11 Tx 13 dRo 15 dDo 16 duRo 17 duDo 18 dPo 19 Nfit 20 chisq 22
     295
     296  foreach f uD uR P
     297    set d$f = $f\o - $f\x
     298  end
     299  set sdR = 3600*(Ro - Rstk)*dcos(Do)
     300  set sdD = 3600*(Do - Dstk)
     301
     302  set Rf = Rx - uRx*(Tx - $Tref)/3600.0/dcos(Dx)
     303  set Df = Dx - uDx*(Tx - $Tref)/3600.0
     304
     305  set dRf = 3600*(Ro - Rf)*dcos(Dx)
     306  set dDf = 3600*(Do - Df)
     307
     308  # echo "fit scatter"
     309  vstat -q dRf;  $value_dRf =  $SIGMA*1000
     310  vstat -q dDf;  $value_dDf =  $SIGMA*1000
     311  vstat -q duR;  $value_duR =  $SIGMA*1000
     312  vstat -q duD;  $value_duD =  $SIGMA*1000
     313  vstat -q dP ;  $value_dP  =  $SIGMA*1000
     314
     315  vstat -q sdR;  $value_sdR =  $SIGMA*1000
     316  vstat -q sdD;  $value_sdD =  $SIGMA*1000
     317
     318  # echo "fit errors"
     319  vstat -q dRo;  $value_dRo  = $MEAN*1000; $sigma_dRo  = $SIGMA*1000
     320  vstat -q dDo;  $value_dDo  = $MEAN*1000; $sigma_dDo  = $SIGMA*1000
     321  vstat -q duRo; $value_duRo = $MEAN*1000; $sigma_duRo = $SIGMA*1000
     322  vstat -q duDo; $value_duDo = $MEAN*1000; $sigma_duDo = $SIGMA*1000
     323  vstat -q dPo ; $value_dPo  = $MEAN*1000; $sigma_dPo  = $SIGMA*1000
     324
     325  # echo "Nfit"
     326  vstat -q Nfit; $mean_Nfit = $MEAN
     327
     328  # echo chisq
     329  vstat -q chisq; $mean_chisq = $MEAN
     330
     331  fprintf "%5.2f %5.2f %5.2f %5.2f %5.2f : %5.2f %5.2f %5.2f %5.2f %5.2f : %5.2f %5.2f %5.2f %5.2f %5.2f : %6.2f %5.2f" $value_dRf $value_dDf $value_duR $value_duD $value_dP $value_dRo $value_dDo $value_duRo $value_duDo $value_dPo $sigma_dRo $sigma_dDo $sigma_duRo $sigma_duDo $sigma_dPo $mean_Nfit $mean_chisq
     332 
     333  delete -q dRbin oRbin
    263334  for i 0 100
    264335    subset tmp = dRf if (Npts == $i)
Note: See TracChangeset for help on using the changeset viewer.