;+ ; NAME: ; vu_insitu_raw ; PURPOSE: ; Makes time series plots of data values at the location of Earth ; from the raw tomography files (in particular from time-dependent tomography) ; CATEGORY: ; sat/idl/util/vupack ; CALLING SEQUENCE: PRO vu_insitu_raw, hdr, ff , $ ut0 = ut0 , $ fixed_time = fixed_time , $ cast_offset = cast_offset , $ ntrail = ntrail , $ raw_time = raw_time , $ last_raw = last_raw , $ silent = silent , $ _extra = _extra ; INPUTS: ; hdr ; ff ; OPTIONAL INPUT PARAMETERS: ; ut0 = ut0 scalar ; type: float ; array[1]; type: time structure ; forecast time on which to center the time series ; specified as Carrington time or time structure in UT ; (passed to vu_insitu) ; cast_offset=cast_offset ; scalar : type: float ; array[1]; type: time difference structure ; fore- or aftcast time relative to the forecast time ; (i.e. the time at which the tomography was run) ; specified in days or as a time difference structure. ; If 'fixe_time' is not specified then ; cast_offset=0.0 is assumed. ; fixed_time=fixe_time ; scalar ; type: float ; array[1]; type: time structure ; Time for which all reconstructed data are requested ; (as Carrington time or time structure) ; raw_time=raw_time scalar ; type: float ; array[1]; type: time structure ; Forecast time for which to plot time series ; (the nearest actual forecast time is used). ; /last_raw specified instead of 'raw_time'. Selects the most recent ; tomographic run. ; OUTPUTS: ; (plot to screen, printer or image file) ; INCLUDE: @compile_opt.pro ; On error, return to caller ; CALLS: ; InitVar, vu_select, vu_get, IsTime, TimeGet, TimeUnit, IsType ; BadValue, vu_gettime, vu_insitu, Carrington ; PROCEDURE: ; The time axis on the plot is the forecast time stored in the headers ; (retrieved with vu_get(hdr, /forecast_time). ; MODIFICATION HISTORY: ; SEP-2002, Paul Hick (UCSD/CASS; pphick@ucsd.edu) ;- InitVar, silent , 0 InitVar, last_raw, /key at_fixed = IsType(fixed_time, /defined) at_raw = IsType(raw_time , /defined) or last_raw uday = TimeUnit(/day) umin = TimeUnit(/minute) ; First get a list of all files. The output rots, file and imark arrays ; are sorted in time and marker values. all_times = vu_select(all_files, _extra=_extra, /nohdr, imark=all_markers, $ count=all_count, silent=silent) IF all_count EQ 0 THEN RETURN ; 'all_markers' will not exist when reading files from corotating tomography corotating = IsType(all_markers, /undefined) CASE corotating OF 0: BEGIN ; Pick up list of unique markers. Read one header for each marker ; value and extract the forecast times. tmp = all_files[uniq(all_markers, sort(all_markers))] forecast_times = vu_select(tmp, imark=forecast_markers, $ count=forecast_count, silent=silent) forecast_times = vu_get(forecast_times, /forecast_time) ; vu_get returns arrays sorted on time coded into file name, i.e. ; forecast_times and forecast_markers are NOT sorted on forecast time ; (same as sorted on marker value). i = sort(forecast_markers) forecast_markers = forecast_markers[i] forecast_times = forecast_times [i] END 1: BEGIN ; For corotating tomography there is only one reconstructed time (same as ; the forecast time). There are no markers so we create dummy arrays. forecast_times = vu_select(all_files, count=forecast_count, silent=silent) forecast_times = vu_get(forecast_times, /forecast_time) all_markers = lindgen(forecast_count) forecast_markers = all_markers END ENDCASE ; forecast_times and forecast_marker are now an ordered list of forecast times ; with its associated marker value. ; 'grid' gets filled with an index into the all_times array ; For files from corotating tomography grid will be a square matrix ; (recon_count = forecast_count) with only the diagonal filled. IF IsType(ntrail, /defined) THEN BEGIN ; If ntrail is specified only retain the last ntrail marker values ; (i.e. the last ntrail tomographic reconstructions) ntrail = ntrail < forecast_count forecast_times = forecast_times [forecast_count-ntrail:forecast_count-1] forecast_markers = forecast_markers[forecast_count-ntrail:forecast_count-1] forecast_count = ntrail tmp = where_common(all_markers, forecast_markers, count=all_count) IF all_count EQ 0 THEN RETURN all_files = all_files [tmp] all_times = all_times [tmp] all_markers = all_markers[tmp] ENDIF ; Get a list of reconstruction times. recon_times = all_times[uniq(TimeOp(/subtract,all_times,all_times[0],uday))] recon_count = n_elements(recon_times) IF (where(indgen(recon_count) NE sort(TimeOp(/subtract,recon_times,recon_times[0],uday))))[0] NE -1 then message, 'oops' grid = make_array(/long, dim=[recon_count, forecast_count], value=-1) marker_offset = forecast_markers[0] all_markers = all_markers-marker_offset loc = lindgen(all_count) FOR i=0,recon_count-1 DO BEGIN ; Loop over all reconstruction times tmp = where(TimeOp(/subtract,all_times,recon_times[i],uday) eq 0) grid[i,all_markers[tmp]] = loc[tmp] ENDFOR IF at_raw THEN BEGIN IF last_raw THEN raw_time = forecast_times[n_elements(forecast_times)-1] tmp = min( abs( TimeOp(/subtract,forecast_times,raw_time,uday) ), i) raw_time = forecast_times [i] raw_marker = forecast_markers[i] good = where(grid[*,i] NE -1, ngood) IF ngood GT 0 THEN BEGIN IF IsType(hdr, /undef) THEN $ hdr = vu_select(all_files[ grid[good,i] ], /read, ff=ff, silent=silent, _extra=_extra) tm_label = '@'+TimeGet(raw_time, /ymd, from=uday,upto=umin)+' UT' InitVar, ut0, raw_time message, /info, 'time series for forecast run '+tm_label+ ' (marker'+strcompress(raw_marker*(1-corotating))+')' vu_insitu, hdr, ff, ut0=ut0, tm_label=tm_label, silent=silent, _extra=_extra ENDIF ENDIF ELSE BEGIN CASE at_fixed OF 0: BEGIN ; Aft/forecast InitVar, cast_offset, 0d0 IF IsTime(cast_offset) THEN cast_offset = TimeGet(cast_offset, /diff, /full, uday, /scalar) ; There are 'forecast_count' forecast times; add cast_offset (positive for forecasts; ; negative for aftcasts). ut = TimeOp(/add,forecast_times, TimeSet(/diff,day=cast_offset)) tm_label = 't'+(['-','+'])[cast_offset GE 0]+strcompress(string(abs(cast_offset),format='(F10.2)'),/rem)+' days' message, /info, 'time series for '+(['aft','fore'])[cast_offset GE 0]+'cast '+tm_label END 1: BEGIN ; Evolution at fixed time ; Replicate fixed time 'fixed_time' ut = replicate(Carrington(fixed_time,/get_time), forecast_count) tm_label = 'fixed @ '+TimeGet(ut[0], /ymd, upto=umin)+ $ ' UT ('+string(Carrington(ut[0]), format='(F9.4)')+')' message, /info, 'time series at '+tm_label END ENDCASE CASE corotating OF 0: BEGIN ; We need to know where times 'ut' fall w.r.t. to the reconstruction ; times. Do this by linear interpolation on an index array. The result ; 'ut_indx' is an array of indices into the rows (1st dim) of 'grid'. ut_indx = TimeInterpol(findgen(recon_count), recon_times, ut) ; Flag as bad anything outside the valid index range for ; 'grid' [0,recon_count-1) tmp = where(ut_indx LT 0 OR ut_indx GT recon_count-1) IF tmp[0] NE -1 THEN ut_indx[tmp] = BadValue(ut_indx) ; Check for each forecast time (each row in 'grid') whether the time ; ut[i] is bracketed by two data files. If not flag as bad. FOR i=0,forecast_count-1 DO $ IF finite(ut_indx[i]) THEN $ IF grid[floor(ut_indx[i]),i] EQ -1 OR grid[ceil(ut_indx[i]),i] EQ -1 THEN $ ut_indx[i] = BadValue(ut_indx) good = where(finite(ut_indx), ngood) ; Check what's left IF ngood GT 1 AND IsType(hdr, /undefined) THEN BEGIN ut = ut [good] ut_indx = ut_indx[good] ; The floor and ceil function could return the same value (if the ; float value ut_indx[i] is equal to an integer), and the same ; file is read twice. p1 = transpose( [ [floor(ut_indx)], [ceil(ut_indx)] ]) p2 = transpose( [ [good ], [good ] ]) hdr = vu_select(all_files[ grid[p1,p2] ], /read, ff=ff, $ /check, /allow_nrad_changes, silent=silent, _extra=_extra) hdr = reform(hdr, 2,ngood, /overwrite) ff = reform(ff , [(size(ff,/dim))[0:3],2,ngood], /overwrite) ; Intepolate between hdr[0,*] and hdr[1,*]. Result goes to hdr[0,*]. FOR i=0,ngood-1 DO BEGIN IF p1[0,i] NE p2[1,i] THEN BEGIN hdr[0,i] = vu_gettime(hdr[*,i], ff[*,*,*,*,*,i], ut=ut[i], ff=ffi) ff[*,*,*,*,0,i] = ffi ENDIF ENDFOR ; Only retain the interpolated data in hdr[0,*] hdr = reform(hdr[0,*]) ff = reform(ff [*,*,*,*,0,*]) ENDIF END 1: BEGIN ngood = recon_count good = lindgen(ngood) IF ngood GT 1 AND IsType(hdr, /undef) THEN $ ; Read files on diagonal of 'grid' array hdr = vu_select(all_files[ grid[good,good] ], /read, ff=ff, $ /check, /allow_nrad_changes, silent=silent, _extra=_extra) END ENDCASE CASE ngood OF 0: 1: message, /info, 'only one data point left; no plot made' ELSE: begin ; For a fore/aftcast time series for time-dependent tomography 'hdr' ; and 'ff' now contain the data arrays at the fore/aftcast times ; (i.e. at the forecast times with the aft/forecast offset ; added/subtracted. Not setting 'ts_all' makes vu_insitu treat these ; data as a regular time-dependent time series. ; For a time-dependent fixed-time time series 'hdr' and 'ff' contain ; the data arrays at the fixed time 'ut' for a number of different ; tomography runs. Setting 'ts_all' (in combination with ; /use_forecast_time) forces vu_insitu to loop over over all ; hdr[i],ff[i] arrays separately and extract data values at Earth at ; times ut[i]. ; For a corotating fore/aftcast time series 'hdr' and 'ff' contain the ; data array at the forecast times, while 'ut' contains the times at ; the forecast times with the fore/aftcast offset added/subtracted. ; Setting 'ts_all' again forces vu_insitu to loop over over all ; hdr[i],ff[i] arrays separately and extract data values at Earth at ; times ut[i]. IF corotating OR at_fixed THEN ts_all = ut ; Set up time at center of plot (best bases on the forecast times, ; since these will be used for the time axis). InitVar, ut0, forecast_times[ forecast_count/2 ] vu_insitu, hdr, ff, ut0=ut0, ts_all=ts_all, /use_forecast_time, $ cast_offset=cast_offset, tm_label=tm_label, _extra=_extra, silent=silent END ENDCASE ENDELSE RETURN & END