IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 39558


Ignore:
Timestamp:
Apr 28, 2016, 11:09:46 AM (10 years ago)
Author:
eugene
Message:

add datan2, set atan2 to return radians as it should

Location:
trunk/Ohana/src/opihi
Files:
1 added
7 edited

Legend:

Unmodified
Added
Removed

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

    r33662 r39558  
    1919  CastVector (window, OPIHI_INT);
    2020
    21   // minimize (flux1 - flux2*alpha) in window defined by mask
     21  // minimize (flux2 - flux1*A - B) in window defined by mask
    2222  // note that the mask is a SELECTION mask not an EXCLUSION mask
    2323
    24   double F12 = 0.0;
    25   double F22 = 0.0;
     24  double R = 0.0, F1 = 0.0, F11 = 0.0, F12 = 0.0, F2 = 0.0;
    2625  for (i = 0; i < flux1->Nelements; i++) {
    27     if ((mask & window->elements.Int[i]) == 0) continue;
     26    // if ((mask & window->elements.Int[i]) == 0) continue;
    2827    double weight = 1.0 / (SQ(dflux1->elements.Flt[i]) + SQ(dflux2->elements.Flt[i]));
     28    R   += weight;
     29    F1  += flux1->elements.Flt[i] * weight;
     30    F2  += flux2->elements.Flt[i] * weight;
     31    F11 += flux1->elements.Flt[i] * flux1->elements.Flt[i] * weight;
    2932    F12 += flux1->elements.Flt[i] * flux2->elements.Flt[i] * weight;
    30     F22 += flux2->elements.Flt[i] * flux2->elements.Flt[i] * weight;
    3133  }
    3234
    33   double Ao = F12 / F22;
    34   double dA = sqrt(1.0 / F22);
     35  int nterm = 2;
     36  double **b = NULL, **c = NULL;
     37  ALLOCATE (b, double *, nterm);
     38  ALLOCATE (c, double *, nterm);
     39  for (i = 0; i < nterm; i++) {
     40    ALLOCATE (c[i], double, nterm);
     41    ALLOCATE (b[i], double, 1);
     42  }
     43  c[0][0] = F11;
     44  c[1][0] = c[0][1] = F1;
     45  c[1][1] = R;
    3546
    36   int Ndof = -1; // 1 parameter fit
     47  b[0][0] = F12;
     48  b[1][0] = F2;
     49
     50  if (!dgaussjordan (c, b, nterm, 1)) {
     51    gprint (GP_ERR, "failed to fit data : ill-conditioned matrix\n");
     52    goto escape;
     53  }
     54
     55  double Ao = b[0][0];
     56  double dA = sqrt(c[0][0]);
     57  double Bo = b[1][0];
     58  double dB = sqrt(c[1][1]);
     59
     60  for (i = 0; i < nterm; i++) {
     61    free (b[i]);
     62    free (c[i]);
     63  }
     64  free (b);
     65  free (c);
     66
     67  int Ndof = -2; // 2 parameter fit
    3768  double chisq = 0.0;
    3869  for (i = 0; i < flux1->Nelements; i++) {
    3970    if ((mask & window->elements.Int[i]) == 0) continue;
    4071    double weight = 1.0 / (SQ(dflux1->elements.Flt[i]) + SQ(dflux2->elements.Flt[i]));
    41     chisq += SQ(flux1->elements.Flt[i] - Ao * flux2->elements.Flt[i]) * weight;
     72    chisq += SQ(flux1->elements.Flt[i] - Ao * flux2->elements.Flt[i] - Bo) * weight;
    4273    Ndof ++;
    4374  }
    … …  
    4576  double chisqNu = chisq / Ndof;
    4677
     78  // fprintf (stderr, "R: %f, F1: %f, F11: %f, F2: %f, F12: %f\n", R, F1, F11, F2, F12);
    4779  // fprintf (stderr, "Ao: %f +/- %f, chisq: %f, chisq_nu : %f for %d dof\n", Ao, dA, chisq, chisqNu, Ndof);
    4880  set_variable ("Ao", Ao);
    4981  set_variable ("dA", dA);
     82  set_variable ("Bo", Bo);
     83  set_variable ("dB", dB);
    5084  set_variable ("Xv", chisqNu);
    5185  set_variable ("Nd", Ndof);

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

    r34088 r39558  
    5353          X = x; Y = y;
    5454          if ((X > -1) && (X < Nx) && (Y > -1) && (Y < Ny)) {
     55            if (!isfinite(*Vin)) continue;
    5556            Vout[X + Y*Nx] += *Vin;
    5657            Sout[X + Y*Nx] ++;

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

    r39225 r39558  
    149149$(SRC)/tvcontour.$(ARCH).o         \
    150150$(SRC)/tvgrid.$(ARCH).o            \
     151$(SRC)/type.$(ARCH).o              \
    151152$(SRC)/uniq.$(ARCH).o              \
    152153$(SRC)/unsign.$(ARCH).o            \

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

    r39225 r39558  
    137137int tvcontour        PROTO((int, char **));
    138138int tvgrid           PROTO((int, char **));
     139int opihi_type             PROTO((int, char **));
    139140int uniq             PROTO((int, char **));
    140141int unsign           PROTO((int, char **));
    … …  
    314315  {1, "tvcontour",    tvcontour,        "send contour to image overlay"},
    315316  {1, "tvgrid",       tvgrid,           "RA/DEC grid on an image"},
     317  {1, "type",         opihi_type,       "get scalar/vector/matrix type information"},
    316318  {1, "ungridify",    ungridify,        "convert image region to vector triplet"},
    317319  {1, "uniq",         uniq,             "create a uniq vector subset from a vector"},

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

    r18732 r39558  
    197197
    198198/* return TRUE / FALSE if name is an existant variable */
     199// XXX is this function better than IsScalar?
     200# if (0)
    199201int get_variable_exists (char *name) {
    200202 
    … …  
    205207  return (FALSE);
    206208}
     209# endif
    207210
    208211float get_variable_default (char *name, float dvalue) {

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

    r38441 r39558  
    7676    if (!strcmp (argv[i], "max"))    { type = ST_BINARY; strcpy (argv[i], "U"); goto gotit; }
    7777    if (!strcmp (argv[i], "min"))    { type = ST_BINARY; strcpy (argv[i], "D"); goto gotit; }
    78     if (!strcmp (argv[i], "atan2"))  { type = ST_BINARY; goto gotit; }
     78    if (!strcmp (argv[i], "atan2"))  { type = ST_BINARY; strcpy (argv[i], "a"); goto gotit; }
     79    if (!strcmp (argv[i], "datan2")) { type = ST_BINARY; strcpy (argv[i], "d"); goto gotit; }
    7980    if (!strcmp (argv[i], ","))      { type = ST_COMMA; goto gotit; }
    8081

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

    r38553 r39558  
    343343    case '^': VV_FUNC(ST_SCALAR_FLT, pow (*M1, *M2));
    344344    case '@': VV_FUNC(ST_SCALAR_FLT, DEG_RAD*atan2 (*M1, *M2));
    345     case 'a': VV_FUNC(ST_SCALAR_FLT, DEG_RAD*atan2 (*M1, *M2));
     345    case 'd': VV_FUNC(ST_SCALAR_FLT, DEG_RAD*atan2 (*M1, *M2));
     346    case 'a': VV_FUNC(ST_SCALAR_FLT,         atan2 (*M1, *M2));
    346347    case 'D': VV_FUNC(ST_SCALAR_INT, MIN (*M1, *M2));
    347348    case 'U': VV_FUNC(ST_SCALAR_INT, MAX (*M1, *M2));
    … …  
    454455    case '^': SV_FUNC(ST_SCALAR_FLT, pow (M1, *M2));
    455456    case '@': SV_FUNC(ST_SCALAR_FLT, DEG_RAD*atan2 (M1, *M2));
    456     case 'a': SV_FUNC(ST_SCALAR_FLT, DEG_RAD*atan2 (M1, *M2));
     457    case 'd': SV_FUNC(ST_SCALAR_FLT, DEG_RAD*atan2 (M1, *M2));
     458    case 'a': SV_FUNC(ST_SCALAR_FLT,         atan2 (M1, *M2));
    457459    case 'D': SV_FUNC(ST_SCALAR_INT, MIN (M1, *M2));
    458460    case 'U': SV_FUNC(ST_SCALAR_INT, MAX (M1, *M2));
    … …  
    561563    case '^': VS_FUNC(ST_SCALAR_FLT, pow (*M1, M2));
    562564    case '@': VS_FUNC(ST_SCALAR_FLT, DEG_RAD*atan2 (*M1, M2));
    563     case 'a': VS_FUNC(ST_SCALAR_FLT, DEG_RAD*atan2 (*M1, M2));
     565    case 'd': VS_FUNC(ST_SCALAR_FLT, DEG_RAD*atan2 (*M1, M2));
     566    case 'a': VS_FUNC(ST_SCALAR_FLT,         atan2 (*M1, M2));
    564567    case 'D': VS_FUNC(ST_SCALAR_INT, MIN (*M1, M2));
    565568    case 'U': VS_FUNC(ST_SCALAR_INT, MAX (*M1, M2));
    … …  
    652655    case '^': MV_FUNC(pow (*M1, *M2));
    653656    case '@': MV_FUNC(DEG_RAD*atan2 (*M1, *M2));
    654     case 'a': MV_FUNC(DEG_RAD*atan2 (*M1, *M2));
     657    case 'd': MV_FUNC(DEG_RAD*atan2 (*M1, *M2));
     658    case 'a': MV_FUNC(        atan2 (*M1, *M2));
    655659    case 'D': MV_FUNC(MIN (*M1, *M2));
    656660    case 'U': MV_FUNC(MAX (*M1, *M2));
    … …  
    745749    case '^': VM_FUNC(pow (*M1, *M2));
    746750    case '@': VM_FUNC(DEG_RAD*atan2 (*M1, *M2));
    747     case 'a': VM_FUNC(DEG_RAD*atan2 (*M1, *M2));
     751    case 'd': VM_FUNC(DEG_RAD*atan2 (*M1, *M2));
     752    case 'a': VM_FUNC(        atan2 (*M1, *M2));
    748753    case 'D': VM_FUNC(MIN (*M1, *M2));
    749754    case 'U': VM_FUNC(MAX (*M1, *M2));
    … …  
    825830    case '^': MM_FUNC(pow (*M1, *M2));
    826831    case '@': MM_FUNC(DEG_RAD*atan2 (*M1, *M2));
    827     case 'a': MM_FUNC(DEG_RAD*atan2 (*M1, *M2));
     832    case 'd': MM_FUNC(DEG_RAD*atan2 (*M1, *M2));
     833    case 'a': MM_FUNC(        atan2 (*M1, *M2));
    828834    case 'D': MM_FUNC(MIN (*M1, *M2));
    829835    case 'U': MM_FUNC(MAX (*M1, *M2));
    … …  
    910916    case '^': MS_FUNC(pow (*M1, M2));
    911917    case '@': MS_FUNC(DEG_RAD*atan2 (*M1, M2));
    912     case 'a': MS_FUNC(DEG_RAD*atan2 (*M1, M2));
     918    case 'd': MS_FUNC(DEG_RAD*atan2 (*M1, M2));
     919    case 'a': MS_FUNC(        atan2 (*M1, M2));
    913920    case 'D': MS_FUNC(MIN (*M1, M2));
    914921    case 'U': MS_FUNC(MAX (*M1, M2));
    … …  
    986993    case '^': SM_FUNC(pow (M1, *M2));
    987994    case '@': SM_FUNC(DEG_RAD*atan2 (M1, *M2));
    988     case 'a': SM_FUNC(DEG_RAD*atan2 (M1, *M2));
     995    case 'd': SM_FUNC(DEG_RAD*atan2 (M1, *M2));
     996    case 'a': SM_FUNC(        atan2 (M1, *M2));
    989997    case 'D': SM_FUNC(MIN (M1, *M2));
    990998    case 'U': SM_FUNC(MAX (M1, *M2));
    … …  
    10681076    case '^': SS_FUNC(ST_SCALAR_FLT, pow (M1, M2));
    10691077    case '@': SS_FUNC(ST_SCALAR_FLT, DEG_RAD*atan2 (M1, M2));
    1070     case 'a': SS_FUNC(ST_SCALAR_FLT, DEG_RAD*atan2 (M1, M2));
     1078    case 'd': SS_FUNC(ST_SCALAR_FLT, DEG_RAD*atan2 (M1, M2));
     1079    case 'a': SS_FUNC(ST_SCALAR_FLT,         atan2 (M1, M2));
    10711080    case 'D': SS_FUNC(ST_SCALAR_INT, MIN (M1, M2));
    10721081    case 'U': SS_FUNC(ST_SCALAR_INT, MAX (M1, M2));
Note: See TracChangeset for help on using the changeset viewer.