Issue hla 255 - Address issues associated with routines under copyrig…

…ht (#408)

* Use released version of ci-watson from conda, rather than master commit via pip. (#399)

* Final adjustments for combining bias and dark images

* Update comments based on review.

* Replaced Numerical Recipes routines and cleaned up obsolete routines

* ACS and WFC3 updates to eliminate copyrighted code.

* Updated text in documenation files

* Removed erroneous comment in Dates/Updates

Co-authored-by: Matt Rendina <[email protected]>

pkg/acs/calacs/Dates

@@ -1,3 +1,6 @@
+ 24-Jun-2020:  CALACS 10.2.4 Removed/modified/updated routines which used copyrighted
+                             code. Updates implemented Oct 2019, but installation delayed 
+                             for higher priority updates.
  29-Apr-2020   CALACS 10.2.3 Modified "dodqi" to accommodate the A-to-D saturation threshold 
                              which is now defined in the CCDTAB. The threshold is no longer 
                              a constant and has different behavior pre- and post-SM4.

pkg/acs/calacs/History

@@ -1,3 +1,8 @@
+### 24-Jun-2020 - MDD -- Version 10.2.4 
+    Removed/modified/updated routines which used copyrighted code.
+    Updates implemented Oct 2019, but installation delayed for higher 
+    priority updates.
+
 ### 29-Apr-2020 - MDD -- Version 10.2.3 
     Modified "dodqi" to accommodate the A-to-D saturation threshold 
     which is now defined in the CCDTAB. The threshold is no longer 

pkg/acs/calacs/Updates

@@ -1,3 +1,9 @@
+Update for version 10.2.4 - 24-Jun-2020 (MDD)
+    pkg/acs/calacs/Dates
+    pkg/acs/calacs/History
+    pkg/acs/calacs/Updates
+    pkg/acs/calacs/acsrej/piksrt.c
+
 Update for version 10.2.3 - 29-Apr-2020 (MDD)
 	pkg/acs/calacs/Dates
 	pkg/acs/calacs/History

pkg/acs/calacs/acsrej/piksrt.c

@@ -1,49 +1,73 @@
+#include <stdio.h>
+# include "trlbuf.h"
 /*
     PIKSRT  -   Sorting by straight insertion. 
-
-    Sorts an array arr[0...n-1] of float type into ascending numerical 
-    order of its .value elements. arr is replaced on output by its sorted 
+    Sorts an array arrayOfFloats[0...n-1] of float type into ascending numerical 
+    order of its value elements. arrayOfFloats is replaced on output by its sorted 
     rearrangement.
 
+    IPIKSRT -   Sorting by straight insertion. 
+    Same algorithm as PIKSRT, but now an arrayOfInts is along for the ride. It is
+    sorted according to the ranking for the arrayOfFloats elements simultaneously.  
+    There is an assumption the two input arrays are the same size.
+
     Order: n^2
+*/
 
-    From "Numerical Recipes - The Art of Scientific Computing",
-    Press, W.H., Teukolsky, S.A., Vetterling, W.T. and Flannery, B.P., 
-    2nd edition, Cambridge University Press, 1995.  
+void piksrt (float arrayOfFloats[], int arrayLength) {
+    float value;
+    int position;
+    trlmessage (MsgText);
 
-*/
+    /* Loop over the full input array */
+    for (int j = 1; j < arrayLength; j++) {
 
-void piksrt (float arr[], int n) {
+        /* Get the value to be sorted and the candidate position */
+        value = arrayOfFloats[j];
+        position = j;
 
-    int         i, j;
-    float       a;
+        /* Is preceding value larger than the new value? 
+           If so, move the preceding value to the higher ranking position...
+        */
+        while (position > 0 && arrayOfFloats[position-1] > value) {
+            arrayOfFloats[position] = arrayOfFloats[position-1];
+            position--;
+        }
 
-    for (j = 1; j < n; j++) {
-        a = arr[j];
-        i = j - 1;
-        while ((i >= 0) && (arr[i] > a)) {
-            arr[i+1] = arr[i];
-            i--;
+        /* ...and put the new value in the lower ranking position. */
+        if (position != j) {
+           arrayOfFloats[position] = value;
         }
-        arr[i+1] = a;
     }
 }
-void ipiksrt (float arr[], int n, int brr[]) {
-
-    int         i, j;
-    float       a;
-    int        b;
-
-    for (j = 1; j < n; j++) {
-        a = arr[j];
-        b = brr[j];
-        i = j - 1;
-        while ((i >= 0) && (arr[i] > a)) {
-            arr[i+1] = arr[i];
-            brr[i+1] = brr[i]; 
-            i--;
+
+void ipiksrt (float arrayOfFloats[], int arrayLength, int arrayOfInts[]) {
+    float fValue;
+    int iValue;
+    int position;
+    trlmessage (MsgText);
+
+    /* Loop over the input array */
+    for (int j = 1; j < arrayLength; j++) {
+
+        /* Get the value to be sorted and the candidate position */
+        fValue = arrayOfFloats[j];
+        iValue = arrayOfInts[j];
+        position = j;
+
+        /* Is preceding value larger than the new value? 
+           If so, move the preceding value to the higher ranking position...
+        */
+        while (position > 0 && arrayOfFloats[position-1] > fValue) {
+            arrayOfFloats[position] = arrayOfFloats[position-1];
+            arrayOfInts[position] = arrayOfInts[position-1];
+            position--;
+        }
+
+        /* ...and put the new value in the lower ranking position. */
+        if (position != j) {
+           arrayOfFloats[position] = fValue;
+           arrayOfInts[position] = iValue;
         }
-        arr[i+1] = a;
-        brr[i+1] = b;
     }
 }

pkg/acs/calacs/include/acsversion.h

@@ -2,8 +2,8 @@
 #define INCL_ACSVERSION_H
 
 /* This string is written to the output primary header as CAL_VER. */
-#define ACS_CAL_VER "10.2.3 (29-Apr-2020)"
-#define ACS_CAL_VER_NUM "10.2.3"
+#define ACS_CAL_VER "10.2.4 (24-Jun-2020)"
+#define ACS_CAL_VER_NUM "10.2.4"
 
 /* name and version number of the CTE correction algorithm */
 #define ACS_GEN1_CTE_NAME "PixelCTE 2012"

pkg/stis/calstis/cs2/piksrt.c

@@ -1,30 +1,35 @@
+#include <stdio.h>
 /*
-  PIKSRT  -  Sorting by straight insertion. 
+    PIKSRT  -   Sorting by straight insertion. 
+    Sorts an array arrayOfFloats[0...n-1] of float type into ascending numerical 
+    order of its value elements. arrayOfFloats is replaced on output by its sorted 
+    rearrangement.
 
-   Sorts an array arr[0...n-1] of float type into ascending numerical 
-   order of its .value elements. arr is replaced on output by its sorted 
-   rearrangement.
-
-   Order: n^2
+    Order: n^2
+*/
 
-   From "Numerical Recipes - The Art of Scientific Computing",
-   Press, W.H., Teukolsky, S.A., Vetterling, W.T. and Flannery, B.P., 
-   2nd edition, Cambridge University Press, 1995.  
+void piksrt (float arrayOfFloats[], int arrayLength) {
+    float value;
+    int position;
 
-*/
+    /* Loop over the full input array */
+    for (int j = 1; j < arrayLength; j++) {
 
-void piksrt (float arr[], int n) {
+        /* Get the value to be sorted and the candidate position */
+        value = arrayOfFloats[j];
+        position = j;
 
-	int         	i, j;
-	float         	a;
+        /* Is preceding value larger than the new value? 
+           If so, move the preceding value to the higher ranking position...
+        */
+        while (position > 0 && arrayOfFloats[position-1] > value) {
+            arrayOfFloats[position] = arrayOfFloats[position-1];
+            position--;
+        }
 
-	for (j = 1; j < n; j++) {
-	    a = arr[j];
-	    i = j - 1;
-	    while ((i >= 0) && (arr[i] > a)) {
-	        arr[i+1] = arr[i];
-	        i--;
-	    }
-	    arr[i+1] = a;
-	}
+        /* ...and put the new value in the lower ranking position. */
+        if (position != j) {
+           arrayOfFloats[position] = value;
+        }
+    }
 }

pkg/stis/calstis/cs6/idtalg/calstis6idt.c

@@ -1745,69 +1745,65 @@ return: the median.
 }
 
 
+#define D_SWAP(a,b) { double temp=(a);(a)=(b);(b)=temp; }
 
-
-# define  SWAP(a,b)  temp=(a);(a)=(b);(b)=temp;
-
-/*
- *  SELECT  -  Selects k-th smallest value in vector.
- *
- *  From "Numerical Recipes - The Art of Scientific Computing",
- *  Press, W.H., Teukolsky, S.A., Vetterling, W.T. and Flannery, B.P.,
- *  2nd edition, Cambridge University Press, 1995.
- *
- *  Vectors are 1-indexed, thus the calling sequence should be
- *  something as:
+/* Algorithm in the public domain (http://blog.beamng.com/a-faster-selection-algorithm/)
+ * with modifications to maintain the HSTCAL interface and added some error checking.
+ * Note: The code needs more than 2 elements to work.
  *
- *  med = select ((long)(n+1)/2, (long)n, temp-1);
+ * This routine maintains the same interface as the original routine.
+ * Arrays are 1-indexed, thus the calling sequence should be
+ * something as:
+ *  med = select ((long)(n+1)/2, (long)n, array-1);
  *
- *  where `temp' is a standard, 0-indexed, 1-dimensional C array.
+ * k: the K-th smallest value in array, range = [1:length]
+ * length: length of the array
+ * array: the array to search
  */
-
-static double Select (unsigned long k, unsigned long n, double *arr) {
-
-        unsigned long  i, ir, j, l, mid;
-        double         a, temp;
-
-        l = 1;
-        ir = n;
-        for (;;) {
-            if (ir <= l+1) {
-                if (ir == l+1 && arr[ir] < arr[l]) {
-                    SWAP(arr[l], arr[ir])
-                }
-                return arr[k];
-            } else {
-                mid = (l+ir) >> 1;
-                SWAP(arr[mid], arr[l+1])
-                if (arr[l] > arr[ir]) {
-                    SWAP(arr[l], arr[ir])
-                }
-                if (arr[l+1] > arr[ir]) {
-                    SWAP(arr[l+1], arr[ir])
-                }
-                if (arr[l] > arr[l+1]) {
-                    SWAP(arr[l], arr[l+1])
-                }
-                i = l+1;
-                j = ir;
-                a = arr[l+1];
-                for (;;) {
-                    do i++; while (arr[i] < a);
-                    do j--; while (arr[j] > a);
-                    if (j < i) break;
-                    SWAP(arr[i], arr[j])
-                }
-                arr[l+1] = arr[j];
-                arr[j]   = a;
-                if (j >= k) ir = j - 1;
-                if (j <= k) l  = i;
-            }
+static double Select (unsigned long k, unsigned long length, double *array) {
+    unsigned long l=0, m=length-1, i=l, j=m;
+    double x;
+
+    /* The original routine did no checking internally, but if its calling
+     * routine returns 0.0, it is an "OUT_OF_MEMORY" error.  This is kludgy,
+     * but it is the only method to return an error in the current context.
+     */
+    if ((length < 2) || (k <= 0) || (k > length)) {
+        //sprintf (MsgText, "Requested value %d is out of range (1 - %d}", k, length);
+        //trlerror (MsgText)
+        return (0.0);
+    }
+
+    /* Adjust k to be for a zero-indexed array */
+    array += 1;
+    k--;
+
+    while (l < m) {
+        if(array[k] < array[i]) D_SWAP(array[i], array[k]);
+        if(array[j] < array[i]) D_SWAP(array[i], array[j]);
+        if(array[j] < array[k]) D_SWAP(array[k], array[j]);
+
+        x=array[k];
+        while (j > k && i < k) {
+            do i++; while (array[i] < x);
+            do j--; while (array[j] > x);
+
+            D_SWAP(array[i], array[j]);
         }
+        i++; j--;
 
+        if (j < k) {
+            while (array[i] < x) i++;
+            l = i; j = m;
+        }
+        if (k < i) {
+            while (x < array[j]) j--;
+            m = j; i = l;
+        }
+    }
+    return array[k];
 }
-# undef  SWAP
-
+# undef  D_SWAP
 
 
 /*
@@ -1854,7 +1850,6 @@ int msize;		i: size of mline array
 }
 
 
-
 /*
    Adds ghost to NUV MAMA image.
 */