;+ ; NAME: ; PLOT_HT ; ; PURPOSE: ; This procedure is used to display height-time curves. It reads in a ; height-time file created by one of the movie programs and generates ; a plot. ; ; CATEGORY: ; MOVIE ; ; CALLING SEQUENCE: ; PLOT_HT ; ; INPUTS: ; ; OPTIONAL INPUT PARAMETERS: ; Filename: If filename is present then it is used immediately ; ; OUTPUTS: ; A plot is generated on the screen, and optionally a print file is ; generated of the form idlplot.psnnn, where nnn is a sequential ; number. ; ; OPTIONAL OUTPUT PARAMETERS: ; ; COMMON BLOCKS: ; com_xplot_ht ; ; SIDE EFFECTS: ; Initiates the XMANAGER if it is not already running. ; ; RESTRICTIONS: ; ; PROCEDURE: ; The various widgets are set up and registered. The user selects the ; height-time file to be processed. The file is read in and the data ; points plotted. The user is then able to fit the data to polynomial ; functions of degree 1,2, or 3. The plot can be printed. The speeds ; can be saved to a file. ; ; MODIFICATION HISTORY: ; Written by: Scott Hawley, NRL Summer Student, June 1996 ; Version 2 RA Howard, NRL, Modified plot calls to use utplot ; 15 Oct 96 RAH, widgetized ; 27 Oct 96 RAH, Corrected overplots of interpolated values ; Set new window number before plotting ; 08 Nov 96 RAH, Corrected situation if called without argument ; 10 Nov 96 RAH, Corrected Acceleration = 2*fit_coeff ; 11 Nov 96 RAH, Added plot of position angles ; 25 Jun 02 NBR, Put in obsolete notice. ; ; @(#)plot_ht.pro 1.7 06/25/02 LASCO IDL LIBRARY ;- pro do_plot ; ; section to plot the height-time curves and overplot of the curve fits ; common com_xplot_ht,wcontrolbase,wcurrhtlabl,filename,htname,fit, $ wprnbase,winnum,fac_rsec,ts,rs,tsint,degree, $ tint,rint,vel,tai,date_obs,time_obs,comment,pa, $ winsave,last_plot !x.title = 'Time (UT)' !y.title = 'Height (R/Rsun)' !p.title = filename+': '+STRMID(UTC2STR(ts(0)),0,10) !p.multi = [0,1,degree] IF (!d.name NE 'PS') AND (winnum GT -1) THEN wset,winnum UTPLOT,ts,rs,psym=2 OUTPLOT,tsint,rint,psym=0 np = n_elements(ts) !x.title='Height (R/Rsun)' IF (degree gt 1) THEN BEGIN velfit=fit(1:degree)*(findgen(degree)+1.0) ; 1st deriv wrt time vel= poly((tai-tai(0))/3600,velfit)*fac_rsec ; interpolated speeds velint = poly(tint,velfit)*fac_rsec ; interpolated speeds maxv=fix(MAX(round(vel/100.))) maxv=100*(maxv+1) IF (maxv lt 500) THEN maxv=500 !y.title = 'Velocity (km/s)' PLOT,rs,vel,psym=2, YRANGE=[0,maxv] OPLOT,rint,velint !x.range='' IF (degree gt 2) THEN BEGIN accfit=velfit(1:degree-1)*(findgen(degree-1)+1.0) ; second derivative !y.title = 'Acceleration (km/s^2)' PLOT,rint,poly(tint,accfit)*fac_rsec/3600. ENDIF ELSE BEGIN accel = 2*fit(degree)*fac_rsec/3.600 str = STRING(accel,format='(f6.2)')+' m/s^2' xyouts,0.5,0.12,'Acceleration = '+str,/normal,charsize=1.3 xyouts,0.5,0.16,'Position Angle = '+STRING(MEDIAN(pa), $ FORMAT='(f5.0)'),/normal,charsize=1.3 ENDELSE ENDIF ELSE BEGIN speed = fit(degree)*fac_rsec str = STRING(speed,format='(f6.1)')+' km/s' xyouts,0.6,0.16,'Velocity = '+str ,/normal,charsize=1.3 xyouts,0.6,0.20,'Position Angle = '+STRING(MEDIAN(pa), $ FORMAT='(f5.0)'),/normal,charsize=1.3 xyouts,0.6,0.22,comment,/normal,charsize=1.3 ENDELSE !x.title='' !y.title='' !p.title='' !P.multi='' IF (!d.name NE 'PS') AND (winsave GT -1) THEN WSET,winsave RETURN END pro ht_print,lpname ; ; prints the current plot onto the printer designated by lpname ; common com_xplot_ht,wcontrolbase,wcurrhtlabl,filename,htname,fit, $ wprnbase,winnum,fac_rsec,ts,rs,tsint,degree, $ tint,rint,vel,tai,date_obs,time_obs,comment,pa, $ winsave,last_plot IF (last_plot NE 'pa') AND (last_plot NE 'ht') THEN RETURN ps_setup,0 !x.range=[0,30] ; impose a constant axis range for comparison CASE last_plot OF 'pa': pa_plot 'ht': do_plot ; do the plot ENDCASE IF (lpname EQ '') OR (STRLOWCASE(lpname) EQ 'lp') $ THEN ps_setup,1 ELSE ps_setup,1,lpname !x.range='' RETURN END pro ht_read ; ; reads in the height-time file, stores the height and time into ; arrays. Time is converted to tai and the CDS time structure ; common com_xplot_ht,wcontrolbase,wcurrhtlabl,filename,htname,fit, $ wprnbase,winnum,fac_rsec,ts,rs,tsint,degree, $ tint,rint,vel,tai,date_obs,time_obs,comment,pa, $ winsave,last_plot maxpts=2000 timestr='' space='' temp=0.0 degree=1 rs=fltarr(maxpts) tai=dblarr(maxpts) pa=fltarr(maxpts) i=0 ; save zeroth element instr='' date_obs='' time_obs='' comment='' ; ; Read in the height-time file ; Test to see if the first character is a #, which means a comment ; The file should not have any blank lines ; Convert the date/time string to TAI time ; OPENR,lu,htname,/get_lun WHILE (not eof(lu)) DO BEGIN READF,lu,instr IF (strmid(instr,0,1) ne '#') THEN BEGIN i=i+1 READS,instr,temp,timestr,p,FORMAT='(F6.2,A20,f6.0)' rs(i)=temp tai(i)=utc2tai(str2utc(timestr)) pa(i)=p ; position angle ENDIF ELSE BEGIN IF (strmid(instr,0,10) eq '#DATE-OBS:') $ THEN date_obs=strmid(instr,11,10)$ ELSE IF (strmid(instr,0,10) eq '#TIME-OBS:') $ THEN time_obs=strmid(instr,11,8) $ ELSE IF (strmid(instr,0,9) eq '#COMMENT:') $ THEN comment=strmid(instr,10,strlen(instr)) ENDELSE ENDWHILE FREE_LUN,lu ; ; Truncate the array to include only valid points ; Then sort by time ; tai= tai(1:i) rs = rs(1:i) pa = pa(1:i) ind=SORT(rs) rs=rs(ind) tai=tai(ind) ts = TAI2UTC(tai) ; convert tai time to CDS time structure !P.MULTI='' RETURN END pro ht_plot ; ; fits the data to a polynomial of the input degree and plots the ; curve over the current plot ; common com_xplot_ht,wcontrolbase,wcurrhtlabl,filename,htname,fit, $ wprnbase,winnum,fac_rsec,ts,rs,tsint,degree, $ tint,rint,vel,tai,date_obs,time_obs,comment,pa, $ winsave,last_plot erase numpts = n_elements(ts) ; ; always do linear fit to get the extrapolated limb time ; fit=poly_fit((tai-tai(0))/3600,rs,1) rmin = 0 ; minimum solar radius to extrapolate to offset = (rmin-fit(0))/fit(1) ; ; now fit to desired degree ; fit=poly_fit((tai-tai(0))/3600,rs,degree) str='Fit = '+string(fit(0))+' + '+string(fit(1))+'x' IF (degree ge 2) THEN FOR i=2,degree DO $ str=str+' + '+string(fit(i))+'x^'+strtrim(string(i),2) print,str ; ; Compute an array of evenly spaced times beginning at the time ; of projected initiation from rmin ; tint is time in hours from the start time ; t0 = ts(0).mjd t1 = ts(numpts-1).mjd+1 tint = findgen(100)/100.*(t1-t0)*24.+offset ; ; rint is an array of interpolated heights at the evenly spaced times ; rint = poly(tint,fit) ; ; tsint is an array of CDS time structures for the evenly spaced times ; tsint = REPLICATE (ts(0),100) tsint.time = long(tint*1000.*3600.)+tsint(0).time do_plot last_plot='ht' RETURN END pro pa_plot ; ; Plots the position angle vs R ; common com_xplot_ht,wcontrolbase,wcurrhtlabl,filename,htname,fit, $ wprnbase,winnum,fac_rsec,ts,rs,tsint,degree, $ tint,rint,vel,tai,date_obs,time_obs,comment,pa, $ winsave,last_plot !y.title = 'Position Angle (Deg)' !x.title = 'Height (R/Rsun)' !p.title = filename+': '+STRMID(UTC2STR(ts(0)),0,10) !p.multi = [0,1,degree] IF (!d.name NE 'PS') AND (winnum GT -1) THEN wset,winnum PLOT,rs,pa,psym=2,/ynozero !x.title='' !y.title='' !p.title='' !P.multi='' IF (!d.name NE 'PS') AND (winsave GT -1) THEN WSET,winsave last_plot='pa' RETURN END pro ht_save,velname ; ; saves the velocity information into a file ; common com_xplot_ht,wcontrolbase,wcurrhtlabl,filename,htname,fit, $ wprnbase,winnum,fac_rsec,ts,rs,tsint,degree, $ tint,rint,vel,tai,date_obs,time_obs,comment,pa, $ winsave,last_plot IF (velname EQ '') THEN RETURN print,' saving velocities to '+velname OPENW,lu,velname,/GET_LUN,/APPEND np = n_elements(rs) ; ; Compute the speed for each of the input times ; IF (degree GT 1) THEN BEGIN velfit=fit(1:degree)*(findgen(degree)+1.0) ; 1st deriv wrt time vel = poly((tai-tai(0))/3600.,velfit)*fac_rsec ENDIF ELSE BEGIN vel = fltarr(np)+fit(1)*fac_rsec ENDELSE FOR i=0,np-1 DO PRINTF,lu,rs(i),vel(i),FORMAT='(F12.2,F12.1)' FREE_LUN,lu RETURN END ;------------------------------------------------------------------------------ ; procedure plot_ht_ev ;------------------------------------------------------------------------------ ; This procedure processes the events being sent by the XManager. ; This is the event handling routine for the plot_ht widget. It is ; responsible for dealing with the widget events such as mouse clicks on ; buttons in the plot_ht widget. The tool menu choice routines are ; already installed. This routine is required for the plot_ht widget to ; work properly with the XManager. ;------------------------------------------------------------------------ PRO plot_ht_ev, event common com_xplot_ht,wcontrolbase,wcurrhtlabl,filename,htname,fit, $ wprnbase,winnum,fac_rsec,ts,rs,tsint,degree, $ tint,rint,vel,tai,date_obs,time_obs,comment,pa, $ winsave,last_plot WIDGET_CONTROL, event.id, GET_UVALUE = eventval ;find the user value ;of the widget where ;the event occured ;print,eventval CASE eventval OF ; here is where you would add the actions for your events. Each widget ; you add should have a unique string for its user value. Here you add ; a case for each of your widgets that return events and take the ; appropriate action. 'prnplt': BEGIN ; Print the plot to the printer field_base = widget_base(wprnbase) lpname = cw_field(field_base,/frame,/return_events, $ /string,title='Enter Printer Name',/column) widget_control,field_base,/realize lpn = widget_event(field_base) ; print,lpn.value widget_control,/destroy,field_base ht_print,lpn.value END 'readht': BEGIN htname = pickfile (title='Select Height-Time File') parts = str_sep(htname,!delimiter) filename = parts(n_elements(parts)-1) ; ; Update the current file name ; widget_control,wcurrhtlabl,set_value=filename ; ; read in the file and plot the h-t values ; ht_read degree = 1 ht_plot END '1st': BEGIN ; fit to a linear function sz = SIZE(filename) IF (sz(1) EQ 0) THEN filename='' degree = 1 ht_plot END '2nd': BEGIN ; fit to a quadratic sz = SIZE(filename) IF (sz(1) EQ 0) THEN filename='' degree = 2 ht_plot END '3rd': BEGIN ; fit to a 3rd degree polynomial sz = SIZE(filename) IF (sz(1) EQ 0) THEN filename='' degree = 3 ht_plot END 'PA': BEGIN sz = SIZE(filename) IF (sz(1) EQ 0) THEN filename='' pa_plot END 'savevl': BEGIN ; save the speeds to a file velfile = pickfile(/write) ; print,velfile ht_save,velfile END 'quit': WIDGET_CONTROL, event.top, /DESTROY ;There is no need to ;"unregister" a widget ;application. The ;XManager will clean ;the dead widget from ;its list. ELSE: MESSAGE, "Event User Value Not Found" ;When an event occurs ;in a widget that has ;no user value in this ;case statement, an ;error message is shown ENDCASE ;print,eventval END ;=========== end of plot_ht event handling routine task ============ ;------------------------------------------------------------------------- ; procedure plot_ht ;------------------------------------------------------------------------- ; This routine creates the widget and registers it with the XManager. ; This is the main routine for the plot_ht widget. It creates the ; widget and then registers it with the XManager which keeps track of the ; currently active widgets. ;------------------------------------------------------------------------- PRO plot_ht, fname, GROUP = GROUP message,'This routine is obsolete--',/info message,'Please use xplot_ht instead.',/info print,'(Will continue in 5 seconds.)' print wait,5 common com_xplot_ht,wcontrolbase,wcurrhtlabl,filename,htname,fit, $ wprnbase,winnum,fac_rsec,ts,rs,tsint,degree, $ tint,rint,vel,tai,date_obs,time_obs,comment,pa, $ winsave,last_plot usefile = 0 IF (N_PARAMS() EQ 1) THEN IF (fname NE '') THEN usefile=1 IF (!D.WINDOW eq -1) THEN winsave=-1 else winsave=!D.WINDOW fac_rsec = 6.96E5/3600. ; factor to convert R/hour to km/sec last_plot = '' ; plot_ht cannot have multiple copies running. IF(XRegistered("plot_ht") NE 0) THEN RETURN ;only one instance of ;the plot_ht widget ;is allowed. If it is ;already managed, do ;nothing and return ; Create the main base. The ROW keyword arranges the widget visually. plot_htbase = WIDGET_BASE(TITLE = "plot_ht",/row) ; To add other routines or remove any of these, remove them both below ; and in the plot_ht_ev routine. plot_htbase = widget_base(/column, $ TITLE=' PLOT HEIGHT-TIME CURVE') wtopbase = widget_base(plot_htbase,/row) wcontrolbase = widget_base(wtopbase,/column,/frame, $ xpad=10,ypad=10,space=35) wreadhtbase = widget_base(wcontrolbase,/column) wreadhtbutton = widget_button(wreadhtbase,uvalue='readht', $ value='Read Height-Time File') curr_ht=' unknown ' wreadhtlabl = widget_label(wreadhtbase,/align_center, $ value='Current HT File:') wcurrhtlabl = widget_label(wreadhtbase,/align_center,value=curr_ht) wcurvtpbase = widget_base(wcontrolbase,/column) wcurvftlabl = widget_label(wcurvtpbase,/align_center, $ value='Select Degree of Polynomial') wcurvftbase = widget_base(wcurvtpbase,/row) wlinrlabl = widget_button(wcurvftbase,uvalue='1st',value='1st') wquadlabl = widget_button(wcurvftbase,uvalue='2nd',value='2nd') wthrdlabl = widget_button(wcurvftbase,uvalue='3rd',value='3rd') wpalabl = widget_button(wcurvftbase,uvalue='PA',value='PA') wprnbase = widget_base(wcontrolbase,/column,xpad=10,ypad=10,space=10) wprnpltbutton = widget_button(wprnbase,uvalue='prnplt', $ value='Print the Plot') wsavevlbutton = widget_button(wcontrolbase,uvalue='savevl', $ value='Save Speeds To File') wquitbutton = widget_button (wcontrolbase,uvalue='quit',value='Quit') wdrawbase = widget_base(wtopbase,/column,/frame,space=20,xpad=10) wdraw = widget_draw(wdrawbase,xsize=512,ysize=512) ; Typically, any widgets you need for your application are created here. ; Create them and use plot_htbase as their base. They will be realized ; (brought into existence) when the following line is executed. WIDGET_CONTROL, plot_htbase, /REALIZE ;create the widgets ;that are defined winnum = !d.window IF (usefile EQ 1) THEN BEGIN filename = fname ff = FINDFILE(filename) sz = SIZE(ff) IF (sz(0) NE 0) THEN BEGIN widget_control,wcurrhtlabl,set_value=filename htname=filename ht_read degree = 1 ht_plot ENDIF ENDIF XManager, "plot_ht", plot_htbase, $ ;register the widgets EVENT_HANDLER = "plot_ht_ev", $ ;with the XManager GROUP_LEADER = GROUP ;and pass through the ;group leader if this ;routine is to be ;called from some group ;leader. IF (winsave GT -1) THEN WSET,winsave END ;================== end of plot_ht main routine ====================