;+ ; NAME: ; vu_insitu_persist ; PURPOSE: ; Makes time series plots of data values at the location of Earth ; from the raw tomography files (in particular from time-dependent tomography). ; All available 'raw' tomography runs in the specified input directory are used. ; CATEGORY: ; sat/idl/util/vupack ; CALLING SEQUENCE: PRO vu_insitu_persist, hdr, ff , $ ut0 = ut0 , $ cast_offset = cast_offset , $ ntrail = ntrail , $ silent = silent , $ printer = printer , $ _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; default: 0.0 ; array[1]; type: time difference structure ; fore- or aftcast time relative to the forecast time ; (i.e. relative to the time at which the tomography was run) ; specified in days or as a time difference structure. ; ntrail=ntrail scalar; type: integer; default: none ; if set, then only the most recent 'ntrail' tomography are ; used ; 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, vu_insitucurve ; EXAMPLE: ; vu_insitu_persist,path=filepath(root=getenv('NAGOYA'),sub=['fast','raw'],'nice'),filter='nv3h*',delt=0.25,source='acesw',/correlation,/timeseries ; PROCEDURE: ; The time axis on the plot is the forecast time stored in the headers ; (retrieved with vu_get(hdr, /forecast_time). ; ; Works similar to the fore/aftcast plots made by vu_insitu_raw, but instead ; of using the tomography forecast at t+cast_offset and plotting it against ; the insitu data at the same times, it plots the tomography difference ; between times t+cast_offset and t against the insitu difference between ; the same times. ; MODIFICATION HISTORY: ; NOV-2011, Paul Hick (UCSD/CASS; pphick@ucsd.edu) ; Based on vu_insitu_raw ;- InitVar, silent , 0 InitVar, printer , /key 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. ; Sorting on marker value is the same as sorting on forecast time. 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 a chronologically ordered list ; of forecast times with associated marker values. ; '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 <= 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 ; Only retain files that are part of the selected ntrail tomography runs 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 IF (where(forecast_markers NE forecast_markers[0]+indgen(forecast_count)))[0] NE -1 THEN BEGIN message, /info, 'forecast runs are not contiguous; one or more marker values are missing' RETURN 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) ; Shift marker values, making 0 the lowest value, so they can be used to index 'grid' ; The implicit assumption is made that forecast_markers is a contiguous set of ; numbers. all_markers -= forecast_markers[0] loc = lindgen(all_count) FOR i=0,recon_count-1 DO BEGIN ; Loop over all reconstruction times ; Find all the tomography files at the reconstruction time recon_times[1]. ; Each tomography run contributes zero or one match. tmp = where(TimeOp(/subtract,all_times,recon_times[i],uday) eq 0) ; Store the location of the matching tomography files. grid[i,all_markers[tmp]] = loc[tmp] ENDFOR InitVar, cast_offset, 1.0d0 IF IsTime(cast_offset) THEN cast_offset = TimeGet(cast_offset, /diff, /full, uday, /scalar) ; There are 'forecast_count' forecast times; ; ut = array[2,forecast_count] with ; ut[0,*] the forecast times ; ut[1,*] the forecast_times with cast_offset added ut = [ $ reform( forecast_times, 1,forecast_count) , $ reform(TimeOp(/add, forecast_times, TimeSet(/diff,day=cast_offset)),1,forecast_count) $ ] 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 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[0,*] LT 0 OR ut_indx[0,*] GT recon_count-1 OR $ ut_indx[1,*] LT 0 OR ut_indx[1,*] 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[0,i]) AND finite(ut_indx[1,i]) THEN $ IF grid[floor(ut_indx[0,i]),i] EQ -1 OR grid[ceil(ut_indx[0,i]),i] EQ -1 OR $ grid[floor(ut_indx[1,i]),i] EQ -1 OR grid[ceil(ut_indx[1,i]),i] EQ -1 THEN $ ut_indx[*,i] = BadValue(ut_indx) good = where(finite(ut_indx[0,*]) AND finite(ut_indx[1,*]), 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). In that case the same ; file is read twice; too bad. ; p1,p2 = array[4,ngood] ; p1[0:1,ngood] bracket the forecast times ; p1[2:3,ngood] bracket the persistence times (forecast+cast_offset) p1 = [ $ floor(ut_indx[0,*]), ceil(ut_indx[0,*]) , $ floor(ut_indx[1,*]), ceil(ut_indx[1,*]) $ ] p2 = transpose( [ [good],[good],[good],[good] ] ) ; p1,p2 = array[2,2*ngood] p1 = reform(p1, 2, 2*ngood, /overwrite) p2 = reform(p2, 2, 2*ngood, /overwrite) ; Read all the files pointed by by grid[p1,p2] ; hdr = array [2,2*ngood] ; ff = array [*,*,*,*,2,2*ngood] hdr = vu_select(all_files[ grid[p1,p2] ], /read, ff=ff, $ /check, /allow_nrad_changes, silent=silent, _extra=_extra, /nosort) hdr = reform(hdr, 2, 2*ngood, /overwrite) ff = reform(ff , [(size(ff,/dim))[0:3],2,2*ngood], /overwrite) ; Intepolate between hdr[0,*] and hdr[1,*]. Result goes to hdr[0,*]. FOR i=0,2*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 (stored in hdr[0,*]) ; hdr = array [ngood] ; ff = array [*,*,*,*,ngood] 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 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 THEN ts_all = ut ; Set up time at center of plot (best based 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, _extra=_extra, silent=silent, /noplot, $ f1=f1, f2=f2, t1=t1, t2=t2, sc_label=sc_label IF n_elements(t1) NE n_elements(t2) OR (where(TimeOp(/subtract,t2,t1,uday) NE 0))[0] NE -1 THEN BEGIN message, /info, 'tomography and insitu times differ (delt NOT set?)' RETURN ENDIF i = 2*indgen(ngood) ; f1[i+1,*] are the tomography data at the fore/aftcast time (forecast time + cast_offset) ; f1[i ,*] are the tomography data at the forecast time ; f2[i+1,*] are the insitu data at the fore/aftcast time (forecast time + cast_offset) ; f2[i ,*] are the insitu data at the forecast time f1 = f1[i+1,*]-f1[i,*] ; Tomography difference f2 = f2[i+1,*]-f2[i,*] ; Insitu difference t1 = t1[i] ; Forecast times t2 = t2[i] ; Forecast times (t2=t1) type = vu_type_insitu('earth_insitu', _extra=_extra, hdr=hdr, type=type) CASE printer OF 0: set_page, /zbuffer, xysize=type[0].xysize, /silent, old_device=old_device, _extra=_extra 1: set_page, /printer, aspect=0.89, old_device=old_device, _extra=_extra ENDCASE CASE n_elements(type) EQ 1 OF 0: data_type = round(total(type.data,2)) 1: data_type = type.data ENDCASE vu_InsituCurve, t1, f1, t2, f2, $ ut0 = utt , $ sc_label = sc_label , $ tm_label = tm_label_ , $ /correlation , $ data_type = data_type , $ fmin = type.Fmin , $ fmax = type.Fmax , $ _extra = _extra , $ plot_type = plot_type , $ ; Output silent = silent display = type[0].display label = type[0].label module = display IF IsType(cast_offset, /defined) THEN $ module += '_'+(['a','f'])[cast_offset ge 0]+string(format='(I3.3)',abs(round(24*cast_offset))) module += label IF NOT plot_type[0] AND plot_type[1] THEN module += '_p' ihdr = forecast_count-1 vu_get_page, hdr[ihdr], module, title='', ut=utt, device=old_device, silent=silent, _extra=_extra END ENDCASE RETURN & END