FUNCTION vu_mean, hdr, ff, weights, tt_weight=tt_weight, userinfo=userinfo, $ hdr_index=hdr_index, last_hdr=last_hdr, marker=marker, ff_mean=ff_mean, $ setbad=setbad, nowrap=nowrap, silent=silent ;+ ; NAME: ; vu_mean ; PURPOSE: ; Calculate a weighted average of a set of tomography arrays ; CATEGORY: ; Tomography ; CALLING SEQUENCE: ; ret_hdr = vu_mean(hdr, ff, weights, tt_weight=tt_weight, ff_mean=ff_mean) ; INPUTS: ; hdr array[n]; type: structure ; header information for all data arrays ; ff array[nLng,nLat,nRad,2,n] ; 3D data arrays (velocity and density, or magnetic field components) ; weights array[n]; type: float ; weights to be applied to data arrays ; tt_weight=tt_weight ; array[1]; type: float or time structure ; time at which mean array is required. ; This time is used to calculated the Carrington range and ; range of interest of the mean array by interpolating on the ; the input values in argument 'hdr'. ; OPTIONAL INPUT PARAMETERS: ; hdr_index=hdr_index ; scalar; type: integer; default: 0 ; 0 <= hdr_index <= n-1; identifies element in 'hdr' on which to base ; the header for the mean data ; /last_hdr if set, the hdr_index = n-1 is used ; /marker ???? maybe useful when processing files for same times from different ; tomography run (the marker values are written into 'userinfo' ; setbad=setbad scalar; type: integer ; if set, then bad values in the matrix nearest to tt_weight ; are also flagged bad in the output matrix. ; /nowrap passed to AlignHeliosphere ; OUTPUTS: ; ret_hdr array[1]; type: structure ; header information for returned mean data. The header is based ; on one of the headers in the input 'hdr' array (determined by ; keywords hdr_index and last_hdr) with some fields updated ; (see PROCEDURE). ; ff_mean=ff_mean array[nLng,nLat,nRad,2]; type: float ; weighted mean data ; if n=1 then the input data are returned unmodified ; OPTIONAL OUTPUT PARAMETERS: ; userinfo=userinfo array[2]; type: string ; contains some descriptive information about the averaging ; process (times of original files and weight factors). ; INCLUDE: @compile_opt.pro ; On error, return to caller ; CALLS: ; InitVar, IsType, BadValue, SuperArray, AlignHeliosphere, vu_set ; TimeGet, TimeUnit, Carrington, vu_get, TimeOp ; RESTRICTIONS: ; The input data need to be arranged chronologically, i.e. ; the hdr.time array needs to be monotic (if href=vu_select= and href=vu_read= are ; used to set up the input 'hdr' and 'ff' arrays this condition is satisfied). ; PROCEDURE: ; > The Carrington range and region of interest for the mean array is calculated by ; linear interpolation on time on the input values in the 'hdr' array. ; ; > The arrays to be combined may refer to different ranges of Carrington variables ; (i.e. with different hdr.array and hdr.roi values). Before calculating the weighted ; mean array a linear interpolation in the longitudinal dimension is used to put all ; arrays on the grid needed for the mean array (this is done by href=AlignHeliosphere=). ; ; > The return header 'ret_hdr' is based on one of the headers in the input 'hdr' array ; (the first one unless /last_hdr or hdr_index is set). The following fields are updated: ; ; ret_hdr.time set to tt_weight ; ret_hdr.array obtained by linear interpolation on hdr.array vs. hdr.time at time tt_weight ; ret_hdr.roi obtained by linear interpolation on hdr.roi vs. hdr.time at time tt_weight ; ret_hdr.file set to blank string (href=vu_write= creates a file name based on hdr_time) ; ret_hdr.time_forecast set to tt_weight ; ret_hdr.time_start set to tt_weight ; ret_hdr.marker if all hdr.marker values are the same that value is retained; ; otherwise set to zero (href=vu_write= appends a non-zero ret_hdr.marker ; value to the file name) ; MODIFICATION HISTORY: ; MAR-2001, Paul Hick (UCSD/CASS) ; JUL-2002, Paul Hick (UCSD/CASS) ; Modified error procedure to avoid divide by zero ; The forecast time in the return header is now also set to tt_weight ; NOV-2002, Paul Hick (UCSD/CASS; pphick@ucsd.edu) ; Added /nowrap keyword ;- InitVar, silent , 0 InitVar, marker , /key InitVar, setbad , /key InitVar, last_hdr, /key time_ww = IsType(tt_weight, /defined) cnt = n_elements(hdr) uday = TimeUnit(/day) IF last_hdr THEN hdr_index = cnt-1 InitVar, hdr_index, 0 ; ret_hdr describes the mean data. It is based on hdr[hdr_index] with a ; number of fields modified. ret_hdr = hdr[hdr_index] IF time_ww THEN BEGIN ut_weight = Carrington(tt_weight,/get_time) ret_hdr = vu_set(ret_hdr, uttime=ut_weight) ENDIF ; If cnt=1 then just return the input data unmodified IF cnt GT 1 THEN BEGIN IF time_ww THEN BEGIN ; Find the array limits and region of interest for the mean data by ; linear interpolation on the values for the input data ; (note that 'interpol' requires the hdr.time array to be monotonic) tt = vu_get(hdr, /uttime) ww = vu_get(hdr, /array) ret_hdr = vu_set(ret_hdr, array=TimeInterpol(ww,tt,ut_weight)) ww = vu_get(hdr, /roi) ret_hdr = vu_set(ret_hdr, roi =TimeInterpol(ww,tt,ut_weight)) ENDIF nLng = vu_get(ret_hdr, /nlng) nLat = vu_get(ret_hdr, /nlat) nRad = vu_get(ret_hdr, /nrad) nDat = (size(ff,/dim))[3] ww = vu_get(ret_hdr, /array) dR = (ww[1]-ww[0])/(nLng-1) ; Fraction of rotation per grid spacing dG = (vu_get(hdr, /array))[0,*]-ww[0] ; Offset between matrices in grid spacings dG = reform(dG) mod 1 dG = dG/dR ; Offset between matrices in grid spacings ff_mean = AlignHeliosphere(ff, -dG, nowrap=nowrap) ;ww = SuperArray(weights, nLng*nLat*nRad*2, /lead) ww = SuperArray(weights, nLng*nLat*nRad*nDat, /lead) ; If tt_weight is specified and /setbad is set, then the matrix nearest ; to tt_weight is used to determine which grid points to flag as bad. IF time_ww AND setbad THEN BEGIN tmp = TimeOp(/subtract, vu_get(hdr, /uttime), ut_weight, uday) tmp = min( abs(tmp), nmin) flag = where( 1-finite(ff_mean[*,*,*,*,nmin]) ) ENDIF ; fww is sum of weights. The weight for bad ff_mean values is assumed to be zero. ; ff_mean is weighted sum of function values. ; Where all values of the summed function values are bad, both fww and ff_mean ; will be zero. fww = total(finite(ff_mean)*ww,5) ff_mean = total(ff_mean*ww,5,/nan) ; The following can be done in one statement: ff_mean = ff_mean/fww, ; but that will generate an error message if fww is zero anywhere. ;ff_mean = ff_mean/fww ww = where(fww NE 0.0) IF ww[0] NE -1 THEN ff_mean[ww] = ff_mean[ww]/fww[ww] ww = where(fww eq 0.0) IF ww[0] NE -1 THEN ff_mean[ww] = BadValue(ff_mean) IF IsType(flag, /defined) THEN $ IF flag[0] NE -1 THEN $ ff_mean[flag] = BadValue(ff_mean) userinfo = 'average of '+strcompress(cnt, /rem)+' files' NCoff = floor((vu_get(ret_hdr, /array))[0]) CASE marker OF 0: userinfo = [userinfo, 'times ' +strjoin(TimeGet(vu_get(hdr, /uttime),/ydoy), ' ')] 1: userinfo = [userinfo, 'markers '+strjoin(strcompress(vu_get(hdr, /marker), /rem), ' ')] ENDCASE IF (where(weights NE 1.0))[0] NE -1 THEN $ userinfo = [userinfo, 'weights: '+strjoin(strcompress(weights, /rem), ' ')] ENDIF ELSE BEGIN ff_mean = ff IF silent LE 1 THEN message, /info, $ 'no mean calculated, single file '+TimeGet(vu_get(hdr, /uttime),/ymd) ENDELSE ; More header modifications ; Clear the file name to avoid accidents with overwriting existing files ret_hdr = vu_set(ret_hdr, file='') ret_hdr = vu_set(ret_hdr, time_step=0.0) IF time_ww THEN BEGIN ; In most cases time_start and time_forecast are meaningless at this point. ; We set them to ut_weight unless the values are the same for all headers averaged. ww = vu_get(hdr,/start_time) IF n_elements(uniq(ww)) GT 1 THEN $ ret_hdr = vu_set(ret_hdr, start_time=ut_weight) ww = vu_get(hdr,/forecast_time) ww = TimeOp(/subtract,ww,ww[0],uday) IF n_elements(uniq(ww)) GT 1 THEN $ ret_hdr = vu_set(ret_hdr, forecast_time=ut_weight) ENDIF IF NOT time_ww OR n_elements(uniq(vu_get(hdr, /marker))) GT 1 THEN ret_hdr = vu_set(ret_hdr, marker=0) RETURN, ret_hdr & END