;+ ; NAME: ; inms_spectra_calculations - performs arithmetic with spectra ; function inms_spectra_calculations, axSpectra, $ ;; (i) an array of 1 or more spectra xOperand, $ ;; (i) the spectra (or number) to modify with add=add, $ ;; (i) add xOperand or number to axSpectra subtract=subtract, $ ;; (i) subtract multiply=multiply, $ ;; (i) multiply divide=divide, $ ;; (i) divide Sigma0_M=Sigma0_M, $ ;; (i) uses a variance of zero for the minuend structure(s) Sigma0_S=Sigma0_S ;; (i) uses a variance of zero for the operand ; ; PURPOSE: ; subtracts one spectra from an array of one or more spectra ; to remove the background ; ; RETURN VALUE; ; an array of spectra the same size as the axSpectra containing ; the result of subtracting xBackground from each record in axSpectra ; ; ARGUMENTS: ; axSpectra the array of minuend spectra ; xOperand the subtrahend spectra ; add, subtract, multiply, divide operation to perform ; ; CHANGES ; 16 Sep 2004 NAB begat coding ; 27 Sep 2004 NAB changed from l1a data structure to spectra structure ; (see inms_get_spectra for struct. def.) ; 28 Sep 2004 NAB added error calculations. added Sigma0_* keywords. ; 10 Oct 2004 NAB replaced 0 with O and O with 0 in confused variable ; definitions ; 3-Jan-2005 DAG enhanced check of spectra structure validity, no longer ; relies on length (subject to change). Instead uses field ; names. ;- ;============================================================================= ;;creates local copy axResult = axSpectra ;;checks inputs on_error, 2 ;;size call of the array of minuend spectra anSizeM = size(axSpectra) nStructures = anSizeM[0] EQ 1?1:anSizeM[2] ;;size call of the subtrahend spectra anSizeS = size(xOperand) If anSizeS[0] EQ 0 then begin bNumber = 1 xOperand = float(xOperand) ;;float of float is still float anSizeS = size(xOperand) ENDIF ELSE bNumber = 0 if inms_validate_spectra_data(axspectra) eq 0 then begin message, 'Supplied spectra structure '$ +'appears to be invalid' endif ;;reports error if second argument invalid If (bNumber EQ 0 AND inms_validate_spectra_data(xOperand) eq 0) or $ (bNumber EQ 1 AND anSizeS[1] NE 4) $ then message, 'You must supply a number or' + $ ' a valid data structure as the second argument' ;;sees what operation is desired and sends error if more or less than one called. bAdd = keyword_set(add) bSubtract = keyword_set(subtract) bMultiply = keyword_set(multiply) bDivide = keyword_set(divide) If (bAdd +bSubtract + bMultiply + bDivide) NE 1 then message, $ 'You must supply exactly one operation keyword' ;;sees if Sigma0_* keywords are set bSigma0_M = keyword_set(Sigma0_M) bSigma0_S = keyword_set(Sigma0_S) ;;number of data points nDataPoints = n_elements(axSpectra[*, 0]) ;;covariance arrays anC1covar = dblarr(nDataPoints, nStructures) anC2covar = dblarr(nDataPoints, nStructures) ;;assorted minuend arrays anCountsC1_M = axSpectra.anC1counts anCountsC2_M = axSpectra.anC2counts anC1Error_M = bSigma0_M?0:axSpectra.anC1error anC2Error_M = bSigma0_M?0:axSpectra.anC2error anC1Mean_M = axSpectra.anC1error^2 anC2Mean_M = axSpectra.anC1error^2 ;;assorted operand arrays (when spectra rather than # is input) if bNumber then begin anCountsC1_O = xOperand anCountsC2_O = xOperand endif else begin anCountsC1_O = xOperand.anc1counts anCountsC2_O = xOperand.anc2counts anC1Error_O = bSigma0_S?0:xOperand.anC1error anC2Error_O = bSigma0_S?0:xOperand.anC2error anC1mean_O = xOperand.anC1error^2 anC2mean_O = xOperand.anC2error^2 ;;finds covariance FOR nI=0, nStructures - 1 DO begin anC1covar[nI] = $ total((anCountsC1_M[*, nI] - anC1Mean_M[*, nI])$ *(anCountsC1_O - anC1mean_O))/nDataPoints anC2covar[nI] = $ total((anCountsC2_M[*, nI] - anC2Mean_M[*, nI])$ *(anCountsC2_O - anC2mean_O))/nDataPoints ENDFOR anC1covar = sqrt(anC1covar)*bSigma0_M*bSigma0_S anC2covar = sqrt(anC2covar)*bSigma0_M*bSigma0_S endelse ;;addition IF bAdd THEN begin anC1error_O = bNumber?0:anC1error_O anC2error_O = bNumber?0:anC2error_O FOR nI=0, nStructures - 1 DO begin axResult[*, nI].anc1counts = anCountsC1_M[*, nI] + anCountsC1_O axResult[*, nI].anc2counts = anCountsC2_M[*, nI] + anCountsC2_O axResult[*, nI].anc1error = $ sqrt(anC1error_M[*, nI]^2 + anC1error_O^2 + 2*anC1Covar[*, nI]) axResult[*, nI].anc2error = $ sqrt(anC1error_M[*, nI]^2 + anC2error_O^2 + 2*anC2Covar[*, nI]) ENDFOR ENDIF ;;subtraction IF bsubtract THEN begin anC1error_O = bNumber?0:anC1error_O anC2error_O = bNumber?0:anC2error_O FOR nI=0, nStructures - 1 DO begin axResult[*, nI].anc1counts = anCountsC1_M[*, nI] - anCountsC1_O axResult[*, nI].anc2counts = anCountsC2_M[*, nI] - anCountsC2_O axResult[*, nI].anc1error = $ sqrt(anC1error_M[*, nI]^2 + anC1error_O^2 + 2*anC1Covar[*, nI]) axResult[*, nI].anc2error = $ sqrt(anC1error_M[*, nI]^2 + anC2error_O^2 + 2*anC2Covar[*, nI]) ENDFOR ENDIF ;;multiplication IF bMultiply THEN begin anC1error_O = bNumber?xOperand:anC1error_O anC2error_O = bNumber?xOperand:anC2error_O IF ~bNumber THEN begin FOR nI=0, nStructures - 1 DO begin axResult[*, nI].anc1counts = anCountsC1_M[*, nI] * anCountsC1_O axResult[*, nI].anc2counts = anCountsC2_M[*, nI] * anCountsC2_O axResult[*, nI].anC1error = $ anC1error_M[*, nI]^2*anC1counts_O^2 + $ anC1error_O^2*anC1counts_M[*, nI]^2 + $ 2*anC1covar[*, nI]^2 axResult[*, nI].anC2error = $ anC2error_M[*, nI]^2*anC2counts_O^2 + $ anC2error_O^2*anC2counts_M[*, nI]^2 + $ 2*anC2covar[*, nI]^2 ENDFOR ENDIF ELSE begin FOR nI=0, nStructures - 1 DO begin axResult[*, nI].anc1counts = anCountsC1_M[*, nI] * xOperand axResult[*, nI].anc2counts = anCountsC2_M[*, nI] * xOperand ENDFOR axResult.anC1error = (anC1error_M*xOperand)^2 axResult.anC2error = (anC2error_M*xOperand)^2 ENDELSE axResult.anC1error = sqrt(axResult.anC1error) axResult.anC2error = sqrt(axResult.anC2error) ENDIF ;;division IF bDivide THEN begin anC1error_O = bNumber?xOperand:anC1error_O anC2error_O = bNumber?xOperand:anC2error_O IF ~bNumber then begin FOR nI=0, nStructures - 1 DO begin axResult[*, nI].anc1counts = anCountsC1_M[*, nI] / anCountsC1_O axResult[*, nI].anc2counts = anCountsC2_M[*, nI] / anCountsC2_O axResult[*, nI].anC1error = $ anC1error_M[*, nI]^2*anC1counts_O^2 + $ anC1error_O^2*anC1counts_M[*, nI]^2 - $ 2*anC1covar[*, nI]^2 axResult[*, nI].anC2error = $ anC2error_M[*, nI]^2*anC2counts_O^2 + $ anC2error_O^2*anC2counts_M[*, nI]^2 - $ 2*anC2covar[*, nI]^2 ENDFOR ENDIF else begin axResult.anC1error = (anC1error_M/xOperand)^2 axResult.anC2error = (anC1error_M/xOperand)^2 ENDELSE FOR nI=0, nStructures - 1 DO begin axResult[*, nI].anc1counts = anCountsC1_M[*, nI] / anCountsC1_O axResult[*, nI].anc2counts = anCountsC2_M[*, nI] / anCountsC2_O ENDFOR axResult.anC1error = sqrt(axResult.anC1error) axResult.anC2error = sqrt(axResult.anC2error) ENDIF return, axResult END