#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <stdio.h>
#include "pslib.h"
#include "stac.h"

psArray *stacErrorImages(psArray *inputs, // Array of input images
                         float gain,    // Gain, in e/ADU
                         float rn       // Read noise, in e
    )
{
    float invGain = 1.0/gain;           // Inverse square root of gain
    rn /= gain;                         // Read noise in ADU
    psArray *errors = psArrayAlloc(inputs->n);

    psTrace("stac.errors", 1, "Calculating error images....\n");

    // Iterate over the input images
    for (int i = 0; i < inputs->n; i++) {
        psTrace("stac.errors",5,"Working on image #%d\n",i);

        psImage *image = inputs->data[i]; // Pull out the image of interest
        int numRows = image->numRows;   // Number of rows
        int numCols = image->numCols;   // Number of columns
        psImage *error = psImageAlloc(numCols, numRows, PS_TYPE_F32); // The error image

        // Iterate over the pixels
        for (int r = 0; r < numRows; r++) {
            for (int c = 0; c < numCols; c++) {
                // We actually calculate the variance
                error->data.F32[r][c] = image->data.F32[r][c]*invGain + rn*rn;
            }
        }

        // Put image onto the array
        errors->data[i] = error;
    }

    return errors;
}
