pro PlotCompareInsitu, t1, f1, t2, f2, $ ut = ut, $ trange = trange, $ forecast= forecast, $ oneplot = oneplot, $ noplot = noplot, $ nocorr = nocorr, $ density = density, $ speed = speed, $ bradial = bradial, $ btang = btang, $ correlation=correlation,$ mag1 = mag1, $ mag2 = mag2, $ ymin = ymin, $ ymax = ymax, $ upto = upto, $ sc_label= sc_label, $ labels = labels, $ label_t = label_t, $ ; same as labels silent = silent, $ _extra = _extra @compile_opt.pro ; On error, return to caller ;+ ; NAME: ; PlotCompareInsitu ; PURPOSE: ; Plots time series for tomography and in situ data; ; and/or correlation plots of tomography vs. in situ data ; CALLING SEQUENCE: ; PlotCompareInsitu, t1, f1, t2, f2 ; INPUTS: ; t1 array[n]; type: time structure ; times of available data points for tomographic data ; f1 array[n,4]; type: float ; tomographic data array with columns for ; density, velocity, radial and tangential ; magnetic field, respectively ; 'bad data' should be indicated by !values.f_nan ; t2 array[n]; type: time structure ; times of available data point for in situ data ; f2 array[n,6]; type: float ; in situ data array with columns for ; density, velocity, radial, tangential, normal ; magnetic field components, and magnetic field ; strength respectively ; 'bad data' should be indicated by !values.f_nan ; OPTIONAL INPUT PARAMETERS: ; ut=ut scalar; type: integer or time structure ; If absent, then array 't1' defines the timeaxis ; If present, the 't1' determines the time range of the ; time axis, but the time axis will be centered on 'ut'. ; If ut is specified as an integer then it is interpreted ; as the difference between local time and UT in hours ; trange=trange ; array[2]; float or time structure ; defines the range of the time axis ; a float array is interpreted as a time range in days ; centered on ut ; /forecast only used if keyword 'ut' is set. ; Adds a vertical dashed line at the center of the time ; series (marking time 'ut'), and adds the word 'FORECAST' ; to the right of the dashed line. ; ; If ut is absent then the time series reflects the time range of the 't1' array. ; ; /mag1 selects radial and tangential magnetic field components ; /mag2 select magnetic field strength ; ; ymin = ymin scalar, array[2]; type: any; scalar is the same [scalar,scalar] ; ymax = ymax scalar, array[2]; type: any; scalar is the same [scalar,scalar] ; sets minimum and maximum on the two timeseries plotted ; ; /noplot suppress graphics output ; ; oneplot = oneplot ; scalar; type: integer; value: 1,2,3,4 ; can be used to plot only one of four plots ; /density alternative for oneplot=1 (overrides 'oneplot' setting) ; /speed alternative for onelplt=2 (overrides 'oneplot' setting) ; ; sc_label scalar; type: string ; label used to indicate origin of in situ data ; (usually this is obtained from href=Instrument= with ; the keyword /label set). ; labels=labels scalar or array; type: string; default: ['corotating', 'time dependent'] ; labels used for the tomography time series ; ; _extra=_extra IDL graphics keyword passed to plot statement (charsize, thick, etc.) ; CALLS: ; InitVar, lsqNormalFit, TimeXAxis, PlotCurve, Carrington, Instrument ; gridgen, IsTime, TimeUnit, TimeGet, TimeSet, TimeOp, TimeLimits ; PROCEDURE: ; The 'tomo' arrays are usually extracted from a tomography output file ; The 'insitu' arrays is set up by a call to InsituTimeseries ; MODIFICATION HISTORY: ; NOV-1999, Paul Hick (UCSD/CASS) ; JUL-2002, Paul Hick (UCSD/CASS) ; Added keyword /forecast (before this was functionality was ; included with keyword 'ut' ; AUG-2002, Paul Hick (UCSD/CASS; pphick@ucsd.edu) ; Added keyword trange. Also modified the calculation of trange if it ; is not specified. ;- InitVar, noplot , /key InitVar, nocorr , /key InitVar, correlation, /key InitVar, forecast , /key InitVar, mag1 , /key InitVar, mag2 , /key InitVar, speed , /key InitVar, density , /key InitVar, bradial , /key InitVar, btang , /key InitVar, silent , /key mag1 = mag1 or bradial or btang case n_elements(ymax) of 0: ymax = [-1,-1] 1: ymax = ymax*[1,1] else: ymax = ymax[0:1] endcase case n_elements(ymin) of 0: ymin = [-1,-1] 1: ymin = ymin*[1,1] else: ymin = ymin[0:1] endcase if speed or bradial then oneplot = 1+2*correlation if density or btang then oneplot = 2+2*correlation InitVar, oneplot, 0 oneplot = oneplot[0] if oneplot ne 0 then begin oneplot = (oneplot < 4) > 1 correlation = oneplot ge 3 endif ; mag1=mag2=0: density & velocity ; mag1=1: radial and tangential magnetic field ; mag2=2: magnetic field strength ytomo = fltarr(n_elements(t1),2,/nozero) ysitu = fltarr(n_elements(t2),2,/nozero) case 1 of mag1: begin ytomo = f1[*,2:3,*] ysitu = f2[*,2:3] end mag2: begin ytomo = f1[*,2:3,*] ytomo[*,0] = sqrt( total(ytomo*ytomo,2, /nan) ) ; Tomography: Br^2+Bt^2 = B^2 ytomo[*,1] = ytomo[*,0] ysitu = f2[*,2:4] ysitu[*,0] = sqrt( total(ysitu*ysitu,2, /nan) ) ; In situ: Br^2+Bt^2 ( Bn is NOT zero) ysitu[*,1] = f2[*,5] ; In situ: B ysitu = ysitu[*,0:1] end else: begin ytomo = f1[*,0:1,*] ysitu = f2[*,0:1] end endcase if oneplot eq 0 then !p.multi = [0,2,2]*(1-nocorr)+[0,2,1]*nocorr x0t = 0.03 ; 0.65 ;0.4 x0c = 0.60 y0t = 0.90 y0c = 0.15 xw = 0.12 yh = 0.075 xs = 0.20 iq1 = 0 iq2 = 1 iq3 = 0 iq4 = 1 case oneplot of 0: 1: iq2 = 0 2: iq1 = 1 3: iq4 = 0 4: iq3 = 1 endcase if iq1 eq iq2 then iq4 = iq3-1 if iq3 eq iq4 then iq2 = iq1-1 xtitle = 'day of year in'+strcompress( (TimeGet(t1[0],TimeUnit(/year)))[0] ) quant = ['velocity (km s!E-1!N)','density (cm!E-3!N)','B!Drad!N (nT)','B!Dtang!N (nT)', $ '!MS!X(B!Drad!N!E2!N+B!Dtang!N!E2!N) (nT)','|B| (nT)'] InitVar, sc_label, Instrument(/label) if n_elements(label_t ) ne 0 then labels = label_t if n_elements(labels) lt 1 then labels = 'corotating' if n_elements(labels) lt 2 then labels = [labels, 'time dependent' ] multi = size(ytomo) if multi[0] eq 2 then multi = 1 else multi = multi[3] colored = !d.name ne !ThePrinter and !d.name ne 'PS' sc_color = !p.color iq_colors = [100, 180] *colored+!p.color*(1-colored) ; Blue and red colors = [22, 7, 36, 83, 127]*colored+!p.color*(1-colored) if not noplot and colored then loadct, 12 ; Load color table only if colors needed for iq=iq1,iq2 do begin ysc = reform(ysitu[*,iq]) ymap = reform(ytomo[*,iq,0]) ytitle = quant[2*mag1+4*mag2+iq] ok_map = where( finite(ymap) ) if ok_map[0] eq -1 then message, /info, 'no '+labels[0]+' tomography data available' ok_sc = where( finite(ysc ) ) if ok_sc [0] eq -1 then message, /info, 'no '+sc_label+' '+ytitle+' data available' yrange = [0,0] if ymin[iq] ne -1 then $ yrange[0] = (mag1+(1-iq)*mag2)*ymin[iq] $ else if ok_map[0] ne -1 then begin yrange[0] = min(ymap[ok_map],/nan) if ok_sc[0] ne -1 then yrange[0] = min([yrange[0],ysc[ok_sc]],/nan) yrange[0] = (mag1+(1-iq)*mag2)*yrange[0] endif if ymax[iq] ne -1 then $ yrange[1] = ymax[iq] $ else if ok_map[0] ne -1 then begin yrange[1] = max(ymap[ok_map],/nan) if ok_sc[0] ne -1 then yrange[1] = max([yrange[1],ysc[ok_sc]],/nan) endif ;if ok_map[0] ne -1 then mapav = mean(ymap[ok_map]) else mapav = !values.f_nan ;if ok_sc [0] ne -1 then scav = mean(ysc [ok_sc ]) else scav = !values.f_nan ;if finite(mapav) or finite(scav) then begin ; title = strcompress(mapav)+strcompress(scav)+strcompress(string(mapav/scav,format='(F10.2)')) ;endif else $ title = '' uday = TimeUnit(/day) case IsType(ut, /defined) of 0: begin title = '' if not IsTime(trange) then trange = TimeLimits(t1, /range) end 1: begin ut = Carrington(ut,/get_time) InitVar, upto, TimeUnit(/min) ;if n_elements(upto) eq 0 then upto = TimeUnit(/min) title = TimeGet(ut,/ymd,upto=upto)+' UT' ; If time range for x-axis is specified as a time structure, use it unmodified. ; If it is a time difference structure, use it with ut. ; If it is defined, but is not a structure, interpret as a time ; time difference in uday, and use it with ut. case 1 of IsTime(trange ): ;IsTime(trange, /diff ): trange = TimeOp(/add, ut, trange) IsType(trange, /defined ): trange = TimeOp(/add, ut, TimeSet(/diff,trange,TimeUnit(/day))) else: begin ; If trange not defined than set it such that ut is centered ; in the plot, and trange brackets all times in t1. trange = TimeOp(/subtract, TimeLimits(t1, /range), ut, TimeUnit(/day)) trange = max(abs(trange))*[-1,1]*( 2*(trange[0] lt 0)-1 ) trange = TimeOp(/add, ut, TimeSet(/diff,trange,TimeUnit(/day))) end endcase end endcase if not noplot then begin TimeXAxis, /exact, trange, uday, /noxtitle, _extra=_extra;, title=title if IsTime(ut) then begin xyouts, mean(!x.crange), !y.crange[1]+.04*(!y.crange[1]-!y.crange[0]), title, align=.5, _extra=_extra if forecast then $ xyouts, total(([0,1]+.25*[1,-1])*!x.crange), total(([0,1]+.5*[1,-1])*!y.crange), 'FORECAST', align=.5, _extra=_extra endif if ok_map[0] ne -1 then begin PlotCurve, /oplot, /newyaxis, t1, ytomo[*,iq,0], ytitle=ytitle,$ yrange=yrange, color=iq_colors[iq], _extra=_extra ; Plot 1st tomography series for i=1,multi-1 do PlotCurve, /oplot, t1, ytomo[*,iq,i], $ color=colors[i-1], _extra=_extra ; Plot other tomography series ;if IsTime(ut) then PlotCurve, /oplot, [ut,ut], !y.crange, linestyle=2 if IsTime(ut) and forecast then begin ; Dashed line down center of plot n = 20 dx = [ut,ut] dy = gridgen(n, range=!y.crange) for i=0,n-1,2 do PlotCurve, /oplot, dx, dy[i:i+1], /silent endif dx = !x.crange[1]-!x.crange[0] dy = !y.crange[1]-!y.crange[0] b = !x.crange[0]+x0t*dx w = xw*dx h = !y.crange[0]+y0t*dy+yh*dy if ok_sc[0] ne -1 then begin ; Overplot s/c data if present PlotCurve, /oplot, t2, ysc, linestyle=2*(1-colored), color=sc_color, _extra=_extra h = h-yh*dy plots , b+[0.,w],h*[1.,1.], linestyle=2*(1-colored), color=sc_color, _extra=_extra xyouts, b+(1+xs)*w,h, sc_label, _extra=_extra endif h = h-yh*dy plots , b+[0.,w],h*[1.,1.], color=iq_colors[iq], _extra=_extra xyouts, b+(1+xs)*w,h, labels[0], _extra=_extra for i=1,multi-1 do begin if total(finite(ytomo[*,iq,i])) ne 0 then begin h = h-yh*dy plots , b+[0.,w],h*[1.,1.], color=colors[i-1], _extra=_extra xyouts, b+(1+xs)*w,h, labels[i], _extra=_extra endif endfor endif endif endfor if not nocorr then begin for iq=iq3,iq4 do begin ymap = ytomo[*,iq,0] ysc = ysitu[*,iq] ok_both = where( finite(ymap) and finite(ysc) ) if ok_both[0] ne -1 then begin range = [min([ysc[ok_both],ymap[ok_both]],/nan),max([ysc[ok_both],ymap[ok_both]],/nan)] title = 'correlation'+strcompress(string(correlate(ymap[ok_both],ysc[ok_both]),format='(F10.3)')) a = lsqNormalFit(ymap[ok_both], ysc[ok_both], $ plot = 1-noplot, $ title = title, $ xtitle = labels[0]+' '+quant[2*mag1+4*mag2+iq], $ ytitle = sc_label +' '+quant[2*mag1+4*mag2+iq], $ xrange = range, $ yrange = range, $ _extra = _extra) if not noplot then oplot, !x.crange,!y.crange, linestyle=2 str1 = 'Fit params: '+quant[2*mag1+4*mag2+iq]+' correlation '+labels[0]+' vs. '+sc_label[0]+ $ ' :'+strcompress(correlate(ymap[ok_both],ysc[ok_both])) str2 = 'Fit params: '+strjoin( strcompress(a[*,0]) ) if not silent then begin print, str1 print, str2 endif ;----- ; Kludge for RedHat 7.0 (cass183). For some reason intercepting the above two ; print statements using a pipe in a bash script doesn't work anymore. ; The following appends the same two lines to the file $TUB/fit_params.txt fit_file = filepath(root=getenv('TUB'),'fit_params.txt') case !version.os_family of 'Windows': begin case !version.release lt 5.4 of 0: spawn, 'echo "'+str1+'" >> '+fit_file+' && echo "'+str2+'" >> '+fit_file, /hide 1: spawn, 'echo "'+str1+'" >> '+fit_file+' && echo "'+str2+'" >> '+fit_file endcase end else: spawn, 'echo "'+str1+'" >> '+fit_file+' && echo "'+str2+'" >> '+fit_file endcase ;------ endif else $ !p.multi[0] = !p.multi[0]-1 endfor endif if oneplot eq 0 then !p.multi = 0 return & end