# include # include # include # include # include # include "msarith.h" /* Library with general utilities for performing arithmetic on NICMOS data groups (look for the twin library n_mathstis.c). There are two basic types of operation: image X image and image X numeric constant. Separate routines are provided for each case. Handling of science and time arrays may depend on pixel content (raw counts or count rate). The (first) input array is modified and holds the output result. The folowing table resumes what is done on each pixel: Operation/ SCI ERR DQ TIME SAMP 2nd operand ADD/IMAGE ADD(*1) COMB.(*1) OR ADD ADD SUB/IMAGE SUB COMB. OR 1st op. 1st op. MULT/IMAGE MULT COMB. OR 1st op. 1st op. DIV/IMAGE DIV COMB. OR 1st op. 1st op. ADD/CONST. ADD COMB. --- --- --- SUB/CONST. SUB COMB. --- --- --- MULT/CONST. MULT COMB. --- MULT(*2) --- DIV/CONST. DIV COMB. --- DIV(*2) --- (*1) count_rate == False -> add raw counts count_rate == True -> translate count rate to counts in both input images; add counts and int. times; translate back to count rate. (*2) count_rate == False -> multiply/divide time array by constant. count_rate == True -> copy time from input image. In the inverse operations ("constant / image" and "constant - image") the DQ, TIME and SAMP arrays are left unchanged. Revision history: --------------- 01 Mar 96 - Implementation (I. Busko) 18 Oct 96 - Set divide-by-zero flag to BADPIX (IB) 11 Nov 96 - STIS support - routines renamed with nn_ prefix (IB) 06 May 97 - Support for inverse division and subtraction (IB) 21 Aug 97 - Propagate divisor's DQ when division-by-zero (IB) */ /* -------------------------------------------------------------------- Error routines specific for this math library. There are two routines, for handling the cases were one or both operands are files. The output DQ is either propagated from the divisor (if the divisor's DQ is non-zero) or set to BADPIX (if the divisor's DQ is zero). ---------------------------------------------------------------------- */ void nn_div0err (SingleNicmosGroup *a, int i, int j, float replace, int verb) { if (verb > 1) { sprintf (ErrText,"Division by zero: pixel %d %d flagged.",i+1,j+1); n_warn (ErrText); } Pix(a->sci.data,i,j) = replace; Pix(a->err.data,i,j) = 0.0F; if (DQPix(a->dq.data,i,j) == 0) DQSetPix(a->dq.data,i,j, BADPIX); } void nn_div0err2 (SingleNicmosGroup *a, SingleNicmosGroup *b, int i, int j, float replace, int verb) { if (verb > 1) { sprintf (ErrText,"Division by zero: pixel %d %d flagged.",i+1,j+1); n_warn (ErrText); } Pix(a->sci.data,i,j) = replace; Pix(a->err.data,i,j) = 0.0F; if (DQPix(b->dq.data,i,j) == 0) DQSetPix(a->dq.data,i,j, BADPIX); else DQSetPix(a->dq.data,i,j, DQPix(b->dq.data,i,j)); } /* -------------------------------------------------------------------- MAIN ROUTINES ---------------------------------------------------------------------- */ /* Add two SingleNicmosGroup quintuplets: (*a) = (*a) + (*b) The science arrays are added together, but the operation can be done in two different ways, depending on the status of the count_rate flag: count_rate = False: it is assumed that the pixels contain raw counts, so they are simply added together. count_rate = True: pixels contain count rate; they are translated back to raw counts, added together and translated again to count rate, with check for zero division. The error arrays are combined following the same procedure above. The sample and time arrays are just added together; the data quality arrays are "or'ed". Division by zero is processed by routine div0err(), and the divZero counter tells the number of pixels where the exception happened. It has to be zeroed on input by the caller. A replace value must be provided in case of division by zero. */ int nn_imageAdd (SingleNicmosGroup *a, SingleNicmosGroup *b, Bool count_rate, int verbose, int *divZero, float replace) { int i, j; /* array indexes */ double aerr; /* a error */ double berr; /* b error */ float time; void nn_div0err (SingleNicmosGroup *, int, int, float, int); errno = 0; for (j=0; j < a->sci.data.ny; j++) { for (i=0; i < a->sci.data.nx; i++) { /* time and sample data are just added together */ time = Pix(a->intg.data,i,j); Pix(a->smpl.data,i,j) = Pix(a->smpl.data,i,j) + Pix(b->smpl.data,i,j); Pix(a->intg.data,i,j) = Pix(a->intg.data,i,j) + Pix(b->intg.data,i,j); /* data quality are or'ed */ DQSetPix(a->dq.data,i,j, DQPix(a->dq.data,i,j) | DQPix(b->dq.data,i,j) ); /* science and error data */ if (count_rate) { /* sci pixels contain count rates */ /* check for division by zero */ if (fabs ((double)time) < DBL_MIN) { nn_div0err (a, i, j, replace, verbose); (*divZero)++; } else { Pix(a->sci.data,i,j) = (Pix(a->sci.data,i,j) * time + Pix(b->sci.data,i,j) * Pix(b->intg.data,i,j)) / Pix(a->intg.data,i,j); aerr = Pix(a->err.data,i,j) * Pix(a->intg.data,i,j); berr = Pix(b->err.data,i,j) * Pix(b->intg.data,i,j); Pix(a->err.data,i,j) = sqrt(aerr*aerr + berr*berr); if (errno) n_sqrterr (&Pix(a->err.data,i,j), verbose); else Pix(a->err.data,i,j) = Pix(a->err.data,i,j) / Pix(a->intg.data,i,j); } } else { /* sci pixels contain raw counts */ Pix(a->sci.data,i,j) = Pix(a->sci.data,i,j) + Pix(b->sci.data,i,j); aerr = Pix(a->err.data,i,j); berr = Pix(b->err.data,i,j); Pix(a->err.data,i,j) = sqrt(aerr*aerr + berr*berr); if (errno) n_sqrterr (&Pix(a->err.data,i,j), verbose); } } } return (0); } /* Add numeric constant to SingleNicmosGroup quintuplet: (*a) = (*a) + constant The constant is added to the science array. The constant error is combined with each pixel of the error array. The remaining arrays are left unchanged. */ int nn_constAdd (SingleNicmosGroup *a, mathConst *constant) { int i, j; /* array indexes */ double aerr; /* a error */ double berr2; /* constant error squared */ errno = 0; berr2 = constant->error * constant->error; for (j=0; j < a->sci.data.ny; j++) { for (i=0; i < a->sci.data.nx; i++) { Pix(a->sci.data,i,j) = Pix(a->sci.data,i,j) + (float)constant->value; aerr = Pix(a->err.data,i,j); Pix(a->err.data,i,j) = sqrt(aerr*aerr + berr2); if (errno) n_sqrterr (&Pix(a->err.data,i,j), 0); } } return (0); } /* Subtract two SingleNicmosGroup quintuplets: (*a) = (*a) - (*b) The science data arrays are subtracted; the error arrays are combined; the data quality arrays are "or'ed". Sample and time arrays are unchanged. */ int nn_imageSub (SingleNicmosGroup *a, SingleNicmosGroup *b) { int i, j; /* array indexes */ double aerr; /* a error */ double berr; /* b error */ errno = 0; for (j=0; j < a->sci.data.ny; j++) { for (i=0; i < a->sci.data.nx; i++) { /* science data */ Pix(a->sci.data,i,j) = Pix(a->sci.data,i,j) - Pix(b->sci.data,i,j); /* error data */ aerr = Pix(a->err.data,i,j); berr = Pix(b->err.data,i,j); Pix(a->err.data,i,j) = sqrt(aerr*aerr + berr*berr); if (errno) n_sqrterr (&Pix(a->err.data,i,j), 0); /* data quality */ DQSetPix(a->dq.data,i,j, DQPix(a->dq.data,i,j) | DQPix(b->dq.data,i,j) ); } } return (0); } /* Subtract numeric constant from SingleNicmosGroup quintuplet: (*a) = (*a) - constant The constant is subtracted from the science array. The constant error is combined with each pixel of the error array. The remaining arrays are left unchanged. */ int nn_constSub (SingleNicmosGroup *a, mathConst *constant) { int i, j; /* array indexes */ double aerr; /* a error */ double berr2; /* constant error squared */ errno = 0; berr2 = (float)(constant->error * constant->error); for (j=0; j < a->sci.data.ny; j++) { for (i=0; i < a->sci.data.nx; i++) { Pix(a->sci.data,i,j) = Pix(a->sci.data,i,j) - (float)constant->value; aerr = Pix(a->err.data,i,j); Pix(a->err.data,i,j) = sqrt(aerr*aerr + berr2); if (errno) n_sqrterr (&Pix(a->err.data,i,j), 0); } } return (0); } /* Subtract SingleNicmosGroup quintuplet from numeric constant: (*a) = constant - (*a) Each pixel from the science array is subtracted from the constant. The constant error is combined with each pixel of the error array. The remaining arrays are left unchanged. */ int nn_constSubInv (SingleNicmosGroup *a, mathConst *constant) { int i, j; /* array indexes */ double aerr; /* a error */ double berr2; /* constant error squared */ errno = 0; berr2 = (float)(constant->error * constant->error); for (j=0; j < a->sci.data.ny; j++) { for (i=0; i < a->sci.data.nx; i++) { Pix(a->sci.data,i,j) = (float)constant->value - Pix(a->sci.data,i,j); aerr = Pix(a->err.data,i,j); Pix(a->err.data,i,j) = sqrt(aerr*aerr + berr2); if (errno) n_sqrterr (&Pix(a->err.data,i,j), 0); } } return (0); } /* Multiply two SingleNicmosGroup quintuplets: (*a) = (*a) * (*b) The science data arrays are multiplied together; the error arrays are combined; the data quality arrays are "or'ed". So far, nothing is done with the time and sample arrays */ int nn_imageMult (SingleNicmosGroup *a, SingleNicmosGroup *b) { int i, j; /* array indexes */ double a_db; /* a value * b error */ double b_da; /* b value * a error */ errno = 0; for (j=0; j < a->sci.data.ny; j++) { for (i=0; i < a->sci.data.nx; i++) { /* error data */ a_db = Pix(a->sci.data,i,j) * Pix(b->err.data,i,j); b_da = Pix(b->sci.data,i,j) * Pix(a->err.data,i,j); Pix(a->err.data,i,j) = sqrt (a_db*a_db + b_da*b_da); if (errno) n_sqrterr (&Pix(a->err.data,i,j), 0); /* science data */ Pix(a->sci.data,i,j) = Pix(a->sci.data,i,j) * Pix(b->sci.data,i,j); /* data quality */ DQSetPix(a->dq.data,i,j, DQPix(a->dq.data,i,j) | DQPix(b->dq.data,i,j) ); } } return (0); } /* Multiply numeric constant into SingleNicmosGroup quintuplet: (*a) = (*a) * constant The science array is multiplied by the constant; the error array is combined pixel-by-pixel with the constant error. Data quality and samples arrays are not modified. The time array may be either multiplied by the constant (if the science pixels are in raw counts) or left unchanged (if the science pixels are in count rate units). */ int nn_constMult (SingleNicmosGroup *a, mathConst *constant, Bool count_rate) { int i, j; /* array indexes */ double a_dct; /* a value * constant error */ double ct_da; /* constant value * a error */ errno = 0; for (j=0; j < a->sci.data.ny; j++) { for (i=0; i < a->sci.data.nx; i++) { /* error */ a_dct = Pix(a->sci.data,i,j) * constant->error; ct_da = constant->value * Pix(a->err.data,i,j); Pix(a->err.data,i,j) = sqrt (a_dct*a_dct + ct_da*ct_da); if (errno) n_sqrterr (&Pix(a->err.data,i,j), 0); /* science */ Pix(a->sci.data,i,j) = Pix(a->sci.data,i,j) * (float)constant->value; /* time */ if (!count_rate) Pix(a->intg.data,i,j) = Pix(a->intg.data,i,j) * (float)constant->value; } } return (0); } /* Divide two SingleNicmosGroup quintuplets: (*a) = (*a) / (*b) The science data arrays are divided; the error arrays are combined; the data quality arrays are "or'ed". So far, nothing is done with the time and sample arrays. Division by zero is processed by routine div0err2(), and the divZero counter tells the number of pixels where the exception happened. It has to be zeroed by the caller. */ int nn_imageDiv (SingleNicmosGroup *a, SingleNicmosGroup *b, int verbose, int *divZero, float replace) { int i, j; /* array indexes */ double asci; /* a science */ double bsci; /* b science */ double bsci2; /* b science squared */ double aerr; /* a error */ double berr; /* b error */ void nn_div0err2 (SingleNicmosGroup *, SingleNicmosGroup *, int, int, float, int); errno = 0; for (j=0; j < a->sci.data.ny; j++) { for (i=0; i < a->sci.data.nx; i++) { asci = Pix(a->sci.data,i,j); bsci = Pix(b->sci.data,i,j); bsci2 = bsci*bsci; if (fabs (bsci2) < DBL_MIN) { /* Division by zero */ nn_div0err2 (a, b, i, j, replace, verbose); (*divZero)++; } else { /* error data */ aerr = Pix(a->err.data,i,j); berr = Pix(b->err.data,i,j); Pix(a->err.data,i,j) = sqrt(aerr*aerr/bsci2 + asci*asci*berr*berr/(bsci2*bsci2)); if (errno) n_sqrterr (&Pix(a->err.data,i,j), verbose); /* science data */ Pix(a->sci.data,i,j) = asci / bsci; /* data quality */ DQSetPix(a->dq.data,i,j, DQPix(a->dq.data,i,j) | DQPix(b->dq.data,i,j) ); } } } return (0); } /* Divide SingleNicmosGroup quintuplet by numeric constant: (*a) = (*a) / constant The science array is divided by the constant; the error array is combined pixel-by-pixel with the constant error. Data quality and samples arrays are not modified. The time array may be either divided by the constant (if the science pixels are in raw counts) or left unchanged (if the science pixels are in count rate units). Division by zero has to be checked by the caller. */ int nn_constDiv (SingleNicmosGroup *a, mathConst *constant, Bool count_rate) { int i, j; /* array indexes */ double asci; /* a science */ double aerr; /* a error */ double cval; /* constant value */ double cval2; /* constant value squared */ double cerr; /* constant error */ double aux; errno = 0; cval = constant->value; cval2 = cval*cval; cerr = constant->error; aux = cerr*cerr / (cval2*cval2); for (j=0; j < a->sci.data.ny; j++) { for (i=0; i < a->sci.data.nx; i++) { /* error */ asci = Pix(a->sci.data,i,j); aerr = Pix(a->err.data,i,j); Pix(a->err.data,i,j) = sqrt(aerr*aerr/cval2 + asci*asci*aux); if (errno) n_sqrterr (&Pix(a->err.data,i,j), 0); /* science */ Pix(a->sci.data,i,j) = asci / cval; /* time */ if (!count_rate) Pix(a->intg.data,i,j) = Pix(a->intg.data,i,j) / (float)constant->value; } } return (0); } /* Divide numeric constant by SingleNicmosGroup quintuplet by: (*a) = constant / (*a) The constant is divided by each pixel in the science array. The error array is combined pixel-by-pixel with the constant error. The remaining arrays are left unchanged. Division by zero is processed by routine div0err(), and the divZero counter tells the number of pixels where the exception happened. It has to be zeroed on input by the caller. A replace value must be provided in case of division by zero. */ int nn_constDivInv (SingleNicmosGroup *a, mathConst *constant, int verbose, int *divZero, float replace) { int i, j; /* array indexes */ double cval; /* constant value */ double asci; /* a science */ double asci2; /* a science squared */ double cerr; /* constant error */ double aerr; /* a error */ void nn_div0err (SingleNicmosGroup *, int, int, float, int); cval = constant->value; cerr = constant->error; errno = 0; for (j=0; j < a->sci.data.ny; j++) { for (i=0; i < a->sci.data.nx; i++) { asci = Pix(a->sci.data,i,j); asci2 = asci*asci; if (fabs (asci2) < DBL_MIN) { /* Division by zero */ nn_div0err (a, i, j, replace, verbose); (*divZero)++; } else { /* error data */ aerr = Pix(a->err.data,i,j); Pix(a->err.data,i,j) = sqrt(cerr*cerr/asci2 + cval*cval*aerr*aerr/(asci2*asci2)); if (errno) n_sqrterr (&Pix(a->err.data,i,j), verbose); /* science data */ Pix(a->sci.data,i,j) = cval / asci; } } } return (0); }