Changeset 39223

Timestamp:

Dec 3, 2015, 4:13:55 PM (11 years ago)

Author:

eugene

Message:

add fitpm_irls, fitplx_irls

Location:

branches/eam_branches/ipp-20151113/Ohana/src/opihi

Files:

• cmd.astro (modified) (1 prop)
• cmd.astro/Makefile (modified) (1 diff)
• cmd.astro/fitplx_irls.c (copied) (copied from trunk/Ohana/src/opihi/cmd.astro/fitplx_irls.c )
• cmd.astro/fitpm_irls.c (modified) (2 diffs)
• cmd.astro/init.c (modified) (2 diffs)
• cmd.data/read_vectors.c (modified) (6 diffs)

branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.astro/Makefile

@@ -35 +35 @@
 $(SRC)/fitplx.$(ARCH).o    \
 $(SRC)/fitpm.$(ARCH).o     \
+$(SRC)/fitpm_irls.$(ARCH).o  \
+$(SRC)/fitplx_irls.$(ARCH).o  \
 $(SRC)/fixwrap.$(ARCH).o           \
 $(SRC)/fixcols.$(ARCH).o           \

branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.astro/fitpm_irls.c

@@ -231 +231 @@
   sigma_ols = 0.0;
   for (i = 0; i < Npts; i++) {
-    rx[i] = X[i] - (T[i] * B[0][0] + B[1][0]);
-    ry[i] = Y[i] - (T[i] * B[2][0] + B[3][0]);
+    rx[i] = X[i] - (T[i] * B[1][0] + B[0][0]);
+    ry[i] = Y[i] - (T[i] * B[3][0] + B[2][0]);
     //    u[i] = r[i] /
     sigma_ols += SQ(rx[i]) + SQ(ry[i]);
@@ -280 +280 @@
     sigma_hat = 0.0;
     for (i = 0; i < Npts; i++) {
-      rx[i] = X[i] - (T[i] * B[0][0] + B[1][0]);
-      ry[i] = Y[i] - (T[i] * B[2][0] + B[3][0]);
+      rx[i] = X[i] - (T[i] * B[1][0] + B[0][0]);
+      ry[i] = Y[i] - (T[i] * B[3][0] + B[2][0]);
       u[i] = sqrt(SQ(rx[i] / dX[i]) + SQ(ry[i] / dY[i]));
     }

branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.astro/init.c

@@ -20 +20 @@
 int fitplx                  PROTO((int, char **));
 int fitpm                   PROTO((int, char **));
+int fitpm_irls              PROTO((int, char **));
+int fitplx_irls             PROTO((int, char **));
 int fixwrap                 PROTO((int, char **));
 int fiximage                PROTO((int, char **));
@@ -87 +89 @@
   {1, "fitplx",      fitplx,       "fit proper motion and parallax"},
   {1, "fitpm",       fitpm,        "fit proper motion only"},
+  {1, "fitpm_irls",  fitpm_irls,   "fit proper motion only using irls method"},
+  {1, "fitplx_irls", fitplx_irls,  "fit proper motion and parallax using irls method"},
   {1, "fixwrap",     fixwrap,      "fix megacam over-wrapped pixels"},
   {1, "fiximage",    fiximage,     "fix pixels in an image by interpolation"},

branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.data/read_vectors.c

@@ -41 +41 @@
   time_t TimeReference;
   int i, j, Nskip, Narg, IsCSV, VERBOSE;
-  int Nbytes, Nstart, NELEM, Nelem, nread;
+  int Nbytes, NELEM, nread;
   char *colstr, *c0, *c1, *extname;
   char varname[1024];  // used as a buffer for the names of string fields
@@ -188 +188 @@
   }
   
+  // we allocate one extra byte into which we never read so there will always be a NULL terminating the string
   ALLOCATE (buffer, char, 0x10001);
   bzero (buffer, 0x10001);
-
-  /* 
-  for (i = 0; i < Nskip; i++) {
-    if (scan_line_maxlen (f, buffer, 0x10000) == EOF) {
-      gprint (GP_ERR, "problem reading file %s\n", filename);
-      read_vectors_cleanup();
-      return FALSE;
-    }
-  }
-  */
 
   int Nline_read = 0; // track number of lines read so far (use to skip lines as well)
@@ -209 +200 @@
   // we treat \n\r pair as a single EOL char to handle mac files:
 
-  Nstart = 0; // location of the last valid byte in the buffer (start filling here)
-  Nelem = 0; // number of valid rows read (vector elements)
-  while (TRUE) {
+  int Nstart = 0; // location of the last valid byte in the buffer (start filling here)
+  int Nelem = 0; // number of valid rows read (vector elements)
+  int EndOfFile = FALSE;
+  while (!EndOfFile) {
     Nbytes = 0x10000 - Nstart;
-    bzero (&buffer[Nstart], Nbytes);
+    // we have allocated one extra byte into which we never read so there will always be a NULL terminating the string
+    bzero (&buffer[Nstart], Nbytes + 1);
     nread = fread (&buffer[Nstart], 1, Nbytes, f);
     if (ferror (f)) {
@@ -219 +212 @@
       break;
     }
-    if (nread == 0) break; // end of the file
-    // nbytes = nread + Nstart;
+    // we still need to parse the rest of the buffer, but there might not be an EOL on the last line
+    if (nread == 0) {
+      EndOfFile = TRUE;
+    }
     
     int bufferStatus = TRUE; 
@@ -227 +222 @@
       c1 = strchr (c0, '\n'); // find the end of this current line (also valid for a Mac: \r\n)
       if (!c1) {
-        c1 = strchr (c0, '\r'); // try \r for non-UNIX files
-      }
-      if (c1 == (char *) NULL) {
+        c1 = strchr (c0, '\r'); // try \r for Windows files
+      }
+      if (!c1) {
         Nstart = strlen (c0);
-        memmove (buffer, c0, Nstart);
-        bufferStatus = FALSE;
-        continue;
+        if (EndOfFile) {
+          // if we have reached EOF, we need to do one last pass in case there is a line without a return
+          c1 = c0 + Nstart;
+          bufferStatus = FALSE;
+        } else {
+          // if we have not reached EOF, we need to shift the buffer to the start of this line and read more data
+          memmove (buffer, c0, Nstart);
+          bufferStatus = FALSE;
+          continue;
+        }
       }
       *c1 = 0; // mark the end of the line 
       Nline_read ++;
 
-      if (Nline_read <= Nskip) { c0 = c1 + 1; continue; }
-
-      if (*c0 == '#') { c0 = c1 + 1; continue; }
-      if (*c0 == '!') { c0 = c1 + 1; continue; }
+      // skip to the next line (but if EOF, do not overrun buffer)
+      if (Nline_read <= Nskip) { if (!EndOfFile) { c0 = c1 + 1; } continue; }
+      if (*c0 == '#')          { if (!EndOfFile) { c0 = c1 + 1; } continue; }
+      if (*c0 == '!')          { if (!EndOfFile) { c0 = c1 + 1; } continue; }
 
       // parse the vectors in this line.  this code is a bit inefficient: each column
@@ -319 +321 @@
         }
       }
-      c0 = c1 + 1;
+      if (!EndOfFile) {
+        c0 = c1 + 1;
+      }
     }
   }