;+
; $Id: animator.pro,v 1.3 2012/02/22 16:51:06 nathan Exp $
;
; Animator is covered in the Sunloop User Guide:
; $SSW/stereo/secchi/idl/display/sunloop/sunloop_user_guide.pdf
;
; written by:   Jeffrey.R.Hall@jpl.nasa.gov
; written for:  Paulett.C.Liewer@jpl.nasa.gov
;               STEREO/SECCHI Project
;
; "Copyright 2012, by the California Institute of Technology. 
; ALL RIGHTS RESERVED. United States Government Sponsorship 
; acknowledged. Any commercial use must be negotiated with the 
; Office of Technology Transfer at the California Institute of 
; Technology.
; This software may be subject to U.S. export control laws. 
; By accepting this software, the user agrees to comply with 
; all applicable U.S. export laws and regulations. User has 
; the responsibility to obtain export licenses, or other 
; export authority as may be required before exporting such 
; information to foreign countries or providing access to 
; foreign persons."
;-

@line3d.pro

FUNCTION get_current_datetime_string
	a=STRSPLIT((SYSTIME()),' ',/EXTRACT)
	CASE a[1] OF
		'Jan' : a[1]='01'
		'Feb' : a[1]='02'
		'Mar' : a[1]='03'
		'Apr' : a[1]='04'
		'May' : a[1]='05'
		'Jun' : a[1]='06'
		'Jul' : a[1]='07'
		'Aug' : a[1]='08'
		'Sep' : a[1]='09'
		'Oct' : a[1]='10'
		'Nov' : a[1]='11'
		'Dec' : a[1]='12'
	ENDCASE
	; Make sure day is 2 digits.
	IF (a[2] LT 10) THEN a[2]='0'+a[2]
	b=STRSPLIT(a[3],':',/EXTRACT)
	RETURN,a[4]+a[1]+a[2]+'_'+b[0]+b[1]+b[2]
END

FUNCTION tick_string, axis, index, value
	;-------------------------------------------------------
	; Plotter tick function to return correct string per tickmark.
	;-------------------------------------------------------
	CALDAT,value,month,day,year,hour,minute,second
	; Fix an odd problem where sometimes time is one minute early.
	minute = CEIL( minute )
	IF ( minute GE 59 ) THEN BEGIN
		minute = 0
		hour = hour + 1
	ENDIF
	; Pad the minute string with a zero, if necessary.
	IF ( minute LT 10 ) THEN minute = '0'+STRTRIM(minute,2) ELSE minute = STRTRIM(minute,2)
	IF ( hour EQ 0 ) THEN BEGIN
		; Date.
		RETURN,STRTRIM(month,2)+'/'+STRTRIM(day,2)+'/'+STRTRIM(year,2)
	ENDIF ELSE BEGIN
		; Time.
		RETURN,STRTRIM(hour,2)+':'+minute
	ENDELSE
END

; example:
; IDL> animator,xyz='Arnaud_sampleXYZ_cloudcme01sep2007_023700.sav.raw.xyz',date='20070901_023700',system='Carrington'

PRO animator, XYZ_FILENAMES=xyz_filenames, DATES=dates, SYSTEM=system

IF KEYWORD_SET(xyz_filenames) THEN import_raw_xyz=1 ELSE import_raw_xyz=0

;-------------------------------------------------------
; Read the data from *.xyz files.  DIALOG_PICKFILE is 
; called for "native" xyz files produced by tiepointer.
; Arbitrary xyz files may be imported on the command line.
;-------------------------------------------------------
IF (import_raw_xyz EQ 0) THEN BEGIN
	xyz_filenames = DIALOG_PICKFILE( /READ, FILTER='*.xyz', /MULTIPLE_FILES )
ENDIF
IF (xyz_filenames[0] EQ '') THEN BEGIN
	; Cancel button was pressed.
	RETURN
ENDIF
IF (STRMID(xyz_filenames[0],STRLEN(xyz_filenames[0])-1,1) EQ PATH_SEP()) THEN BEGIN
	; OK button was clicked without first selecting a file.
	xyz_filenames[0]=''
	RETURN
ENDIF
nfeatures = N_ELEMENTS( xyz_filenames )
xyz_ptrarr = PTRARR( nfeatures )
IF (import_raw_xyz EQ 0) THEN BEGIN
	lat_ptrarr = PTRARR( nfeatures )
	lon_ptrarr = PTRARR( nfeatures )
	rsun_ptrarr = PTRARR( nfeatures )
	system = STRARR( nfeatures )
	telescope = STRARR( nfeatures, 2 )
ENDIF
date = STRARR( nfeatures )
max_count = 0ULL
FOR file = 0ULL, nfeatures - 1 DO BEGIN
	IF (import_raw_xyz EQ 1) THEN BEGIN
		;-------------------------------------------------------
		; Read raw XYZ data, not from tiepioter tool.
		;-------------------------------------------------------
		ascii = READ_ASCII( xyz_filenames[file], COUNT=max_count )
		xyz_ptrarr[file] = PTR_NEW( ascii.(0)[0:2,0:max_count-1] )
		date[file] = dates[file]
	ENDIF ELSE BEGIN
		;-------------------------------------------------------
		; Read embelished XYZ data from tieponter tool.
		;-------------------------------------------------------
		ascii = READ_ASCII( xyz_filenames[file], DATA_START=12 )
		metadata_start = (WHERE( (ascii.(0))[0,*] EQ -9999 ))[0]
		xyz_ptrarr[file] = PTR_NEW( ascii.(0)[4:6,0:metadata_start-1] )
		lat_ptrarr[file] = PTR_NEW( ascii.(0)[7,0:metadata_start-1] )
		lon_ptrarr[file] = PTR_NEW( ascii.(0)[8,0:metadata_start-1] )
		rsun_ptrarr[file] = PTR_NEW( ascii.(0)[9,0:metadata_start-1] )
		max_count = max_count > N_ELEMENTS( (*xyz_ptrarr[file])[0,*] )

		; Get coordinate system, Carrington or Stonyhurst.
		ascii = ''
		OPENR, lun, xyz_filenames[file], /GET_LUN
		FOR line = 0, 12+metadata_start+1 DO BEGIN
			READF, lun, ascii
			IF (line EQ 12+metadata_start+1) THEN BEGIN
				system[file] = ascii
			ENDIF
		ENDFOR
		FREE_LUN, lun

		; Get telescope from left image filename, which is contained inside the xyz file.
		ascii = ''
		OPENR, lun, xyz_filenames[file], /GET_LUN
		FOR line = 0, 12+metadata_start+3 DO BEGIN
			READF, lun, ascii
			; Get wavelength (this only matters if telescope is EUVI.)
			IF (line EQ 11+metadata_start+3) THEN telescope[file,1] = ascii
			; Get telescope.
			IF (line EQ 12+metadata_start+3) THEN BEGIN
				left_filename = FILE_BASENAME((STRSPLIT(ascii,' ',/EXTRACT))[0])
				telescope[file,0] = STRMID((STRSPLIT(left_filename,'_',/EXTRACT))[2],2,2)
				date[file] = STRMID( left_filename, 0, 15 )
			ENDIF
		ENDFOR
		FREE_LUN, lun
	ENDELSE
ENDFOR

IF (import_raw_xyz EQ 0) THEN BEGIN
	IF ((WHERE(system NE system[0]))[0] NE (-1)) THEN BEGIN
		answer=DIALOG_MESSAGE(['WARNING:','NOT ALL XYZ FILES HAVE SAME COORDINATE SYSTEM.','','CONTINUE?'],/QUESTION)
		IF (STRUPCASE(answer) EQ 'NO') THEN BEGIN
			RETURN
		ENDIF
	ENDIF
ENDIF

;-------------------------------------------------------
; Format the data.
; xyz is has 3 dimensions, where the first dimension 
; indexes x, y and z, the second dimension indexes each 
; point, and the third dimension indexes the feature 
; number (where there are multiple features.)
;-------------------------------------------------------
xyz = FLTARR( 3, max_count, nfeatures )
time = DBLARR( nfeatures )
height = DBLARR( nfeatures )
x_count = 0ULL
x_total = 0
y_count = 0ULL
y_total = 0
z_count = 0ULL
z_total = 0
x_count_rmax = 0ULL
x_total_rmax = 0
y_count_rmax = 0ULL
y_total_rmax = 0
z_count_rmax = 0ULL
z_total_rmax = 0
all_x_rmax = DBLARR( nfeatures)
all_y_rmax = DBLARR( nfeatures)
all_z_rmax = DBLARR( nfeatures)
solar_radius_km = 6.96e05

FOR feature = 0, nfeatures - 1 DO BEGIN

	; xyz.
	npoints=N_ELEMENTS((*xyz_ptrarr[feature])[0,*])
	xyz[*,0:npoints-1,feature] = *xyz_ptrarr[feature]

	; remaining variables used only for H-T plots and H-T txt file, not sent to line3d.)
	anytime = ANYTIM2JD( date[feature] )
	time[feature] = anytime.(0) + anytime.(1)
	; height is assumed to be the point farthest from the center of the sun.
	height[feature] = MAX( SQRT( xyz[0,*,feature] ^ 2 + xyz[1,*,feature] ^ 2 + xyz[2,*,feature] ^ 2 ) )
	
	IF (import_raw_xyz EQ 0) THEN BEGIN
		; Convert from degrees to radians.
		lat = (*lat_ptrarr[feature])[*] * CSPICE_RPD()
		lon = (*lon_ptrarr[feature])[*] * CSPICE_RPD()
		; Convert from R/Rsun to Kilometers.
		rsun = (*rsun_ptrarr[feature])[*]
		CSPICE_LATREC, rsun, lon, lat, xxyyzz
		x = REFORM( xxyyzz[0,*] )
		y = REFORM( xxyyzz[1,*] )
		z = REFORM( xxyyzz[2,*] )
		x_count = x_count + N_ELEMENTS( x )
		x_total = TOTAL( [ x_total, x ], /DOUBLE )
		y_count = y_count + N_ELEMENTS( y )
		y_total = TOTAL( [ y_total, y ], /DOUBLE )
		z_count = z_count + N_ELEMENTS( z )
		z_total = TOTAL( [ z_total, z ], /DOUBLE )
        	rmax_ndx = (WHERE(rsun EQ MAX(rsun)))[0]
        	x_count_rmax = x_count_rmax + 1
        	x_total_rmax = TOTAL( [ x_total_rmax, x[rmax_ndx] ], /DOUBLE )
        	y_count_rmax = y_count_rmax + 1
        	y_total_rmax = TOTAL( [ y_total_rmax, y[rmax_ndx] ], /DOUBLE )
        	z_count_rmax = z_count_rmax + 1
	        z_total_rmax = TOTAL( [ z_total_rmax, z[rmax_ndx] ], /DOUBLE )
		; Save these for calculating Standard Deviation.
		all_x_rmax[feature]=x[rmax_ndx]
		all_y_rmax[feature]=y[rmax_ndx]
		all_z_rmax[feature]=z[rmax_ndx]
	ENDIF

ENDFOR

IF (import_raw_xyz EQ 0) THEN BEGIN
	x_mean = x_total / x_count
	y_mean = y_total / y_count
	z_mean = z_total / z_count
	x_mean_rmax = x_total_rmax / x_count_rmax
	y_mean_rmax = y_total_rmax / y_count_rmax
	z_mean_rmax = z_total_rmax / z_count_rmax
	CSPICE_RECLAT, [x_mean,y_mean,z_mean], height_mean, lon_mean, lat_mean
	CSPICE_RECLAT, [x_mean_rmax,y_mean_rmax,z_mean_rmax], height_mean_rmax, lon_mean_rmax, lat_mean_rmax
	lon_mean = lon_mean * CSPICE_DPR()
	lat_mean = lat_mean * CSPICE_DPR()
	lon_mean_rmax = lon_mean_rmax * CSPICE_DPR()
	lat_mean_rmax = lat_mean_rmax * CSPICE_DPR()

	; Calculate rmax Standard Deviation.
	all_lon_rmax = DBLARR(nfeatures)
	all_lat_rmax = DBLARR(nfeatures)
	FOR feature = 0, nfeatures - 1 DO BEGIN
		CSPICE_RECLAT, [ all_x_rmax[feature], all_y_rmax[feature], all_z_rmax[feature] ], height_rmax, lon_rmax, lat_rmax
		all_lon_rmax[feature] = lon_rmax * CSPICE_DPR()
		all_lat_rmax[feature] = lat_rmax * CSPICE_DPR()
	ENDFOR
	lon_sd_rmax = SQRT( ( TOTAL( (all_lon_rmax - MEAN(all_lon_rmax)) ^ 2.0d, /DOUBLE ) ) / nfeatures )
	lat_sd_rmax = SQRT( ( TOTAL( (all_lat_rmax - MEAN(all_lat_rmax)) ^ 2.0d, /DOUBLE ) ) / nfeatures )

	; Free memory.
	PTR_FREE,lat_ptrarr
	PTR_FREE,lon_ptrarr
ENDIF

; Free memory.
PTR_FREE,xyz_ptrarr

;-------------------------------------------------------
; Fix ordering.  This is necessary when XYZ filenames 
; are not time ordered.
;-------------------------------------------------------
sort_ndx = SORT( time )
date = date[sort_ndx]
time = time[sort_ndx]
height = height[sort_ndx]
IF (import_raw_xyz EQ 0) THEN BEGIN
	system = system[sort_ndx]
	telescope = telescope[sort_ndx,*]
ENDIF
xyz_filenames = xyz_filenames[sort_ndx]
xyz = xyz[*,*,sort_ndx]

		check_for_euvi_ndx = WHERE( telescope[*,0] EQ 'eu', euvi_count )
		IF (euvi_count GT 0) THEN BEGIN
			; Append wavelength to eu.  Do this before sort so that sorting includes wavelength.
			telescope[check_for_euvi_ndx,0] = telescope[check_for_euvi_ndx,0] + ' ' + telescope[check_for_euvi_ndx,1]
		ENDIF

IF (nfeatures GT 1) THEN BEGIN
	;-------------------------------------------------------
	; Height-Time plot.
	;-------------------------------------------------------
	km_sun_radius = 695000L
	height = height / km_sun_radius
	xmin=ULONG64( MIN( time ) * 24) / 24.0d
	xmax=ULONG64( MAX( time ) * 24 + 1 ) / 24.0d
	; Reformat date&time string for plot title.
	datetime =	      STRTRIM( FIX( STRMID( date[0], 4, 2 ) ), 2 ) + $	; FIX is used to strip any leading zero from month.
				'/' + STRTRIM( FIX( STRMID( date[0], 6, 2 ) ), 2 ) + $	; FIX is used to strip any leading zero from day.
				'/' + STRMID( date[0], 0, 4 ) + $
				' ' + STRTRIM( FIX( STRMID( date[0], 9, 2 ) ), 2 ) + $	; FIX is used to strip any leading zero from hour.
				':' + STRMID( date[0], 11, 2 )
	; Initialize number of tick intervals to one per hour.
	xticks = ROUND( (xmax - xmin) * 24 + 0 )
	; Restrict the number of ticks intervals to 8 or less.
	interval_length = 1	;hour
	WHILE ( xticks GT 8 ) DO BEGIN
		xticks = xticks / 2
		interval_length = interval_length * 2
	ENDWHILE
	curwin=!D.WINDOW
	WINDOW, /FREE
	win = !D.WINDOW
	PLOT, time, height, TITLE='CME starting at '+datetime, YTITLE='R/Rsun', $
		XTICKFORMAT='tick_string', XRANGE=[xmin,xmax], XTICKS=xticks, XSTYLE=1, $
		COLOR = 0, BACKGROUND = !D.N_COLORS - 1

	IF (import_raw_xyz EQ 0) THEN BEGIN
		;-------------------------------------------------------
		; Overplot symbols on Height-Time plot.
		;-------------------------------------------------------
		; Eliminate second dimension.
		telescope = REFORM( telescope[*,0] )
		tele_sort_ndx = SORT( telescope )
		check_for_euvi_ndx = WHERE( STRMID(telescope[tele_sort_ndx,0],0,2) EQ 'eu', euvi_count )
		check_for_other_ndx = WHERE( STRMID(telescope[tele_sort_ndx],0,2) NE 'eu', other_count )
		IF (euvi_count GT 0) AND (other_count GT 0) THEN BEGIN
			; Push EUVI indices to front.  Otherwise this list is alphabetical [c1,c2,eu,h1,h2].
			tele_sort_ndx = [ tele_sort_ndx[check_for_euvi_ndx], tele_sort_ndx[check_for_other_ndx] ]
		ENDIF
		telescope_sorted = telescope[tele_sort_ndx]
		scope = UNIQ( telescope_sorted )
		time_sorted = time[tele_sort_ndx]
		height_sorted = height[tele_sort_ndx]
		; Prioritize the six useful symbols provided by IDL.
		psyms = [5,6,4,2,7,1]
		FOR tele = 0, N_ELEMENTS( scope ) - 1 DO BEGIN
			ndx = WHERE( telescope_sorted EQ telescope_sorted[scope[tele]] )
			IF (tele LE 5) THEN BEGIN
				OPLOT, time_sorted[ndx], height_sorted[ndx], PSYM=psyms[tele], COLOR = 0
			ENDIF ELSE BEGIN
				; Create more symbols.
				a = FINDGEN(17) * (!PI*2/16.)
				IF (tele EQ 6) THEN USERSYM, COS(a), SIN(a)									; circle
				IF (tele EQ 7) THEN USERSYM, 2*[.5,.85,0,1,.15,.5]-1,2*[1,0,.6,.6,0,1]-1	; 5-point star
				OPLOT, time_sorted[ndx], height_sorted[ndx], PSYM=8, COLOR = 0
			ENDELSE
		ENDFOR
		;-------------------------------------------------------
		; Overplot least squares line on Height-Time plot.
		;-------------------------------------------------------
;		; Manually calculate least squares fit.
;		m = nfeatures
;		x = time - time[0]
;		y = height
;		sumx = TOTAL( x, /DOUBLE )
;		sumy = TOTAL( y, /DOUBLE )
;		sumxsq = TOTAL( x ^ 2, /DOUBLE )
;		sumxy = TOTAL( x * y, /DOUBLE )
;		d = ( m * sumxsq ) - ( sumx ^ 2 )
;		offset = (( sumxsq * sumy ) - ( sumx * sumxy )) / d		; alpha
;		slope = (( m * sumxy ) - ( sumx * sumy )) / d			; beta
		; Use IDL routine to calculate least squares fit.
		slope = (REGRESS( (time - time[0]), height, CONST=offset, /DOUBLE ))[0]
		OPLOT, [time[0],time[nfeatures-1]], [offset,(time[nfeatures-1]-time[0])*slope+offset], $
			LINESTYLE = 1, COLOR = 0
	ENDIF
tzero=time[0]-offset/slope
CALDAT,tzero,m,d,y,h,min,sec
IF (m lt 10) THEN m='0'+STRTRIM(m,2) ELSE m=STRTRIM(m,2)
IF (d lt 10) THEN d='0'+STRTRIM(d,2) ELSE d=STRTRIM(d,2)
IF (h lt 10) THEN h='0'+STRTRIM(h,2) ELSE h=STRTRIM(h,2)
IF (min lt 10) THEN min='0'+STRTRIM(min,2) ELSE min=STRTRIM(min,2)
sec=FIX(sec)
IF (sec lt 10) THEN sec='0'+STRTRIM(sec,2) ELSE sec=STRTRIM(sec,2)
tzero_datetime=STRTRIM(y,2)+m+d+'_'+h+min+sec
	XYOUTS, [.150], [.85], 't_lin_rs = '+tzero_datetime, SIZE = 1.5, COLOR = 0, /NORMAL
	;-------------------------------------------------------
	; Create symbol key for Height-Time plot.
	;-------------------------------------------------------
	PLOT, [0,1], [0,1], /NODATA, /NOERASE, XMARGIN=[0,0], YMARGIN=[0,0], XSTYLE=4, YSTYLE=4
	FOR tele = 0, N_ELEMENTS( scope ) - 1 DO BEGIN
		IF (tele LE 5) THEN BEGIN
			OPLOT, [.150], [.70 - tele / 20.], PSYM = psyms[tele], COLOR = 0
		ENDIF ELSE BEGIN
			; Create more symbols.
			a = FINDGEN(17) * (!PI*2/16.)
			IF (tele EQ 6) THEN USERSYM, COS(a), SIN(a)									; circle
			IF (tele EQ 7) THEN USERSYM, 2*[.5,.85,0,1,.15,.5]-1,2*[1,0,.6,.6,0,1]-1	; 5-point star
			OPLOT, [.150], [.70 - tele / 20.], PSYM = 8, COLOR = 0
		ENDELSE
		CASE telescope_sorted[scope[tele]] OF
			'eu 171'	: tel = 'EUVI 171'
			'eu 195'	: tel = 'EUVI 195'
			'eu 284'	: tel = 'EUVI 284'
			'eu 304'	: tel = 'EUVI 304'
			'c1'		: tel = 'COR1'
			'c2'		: tel = 'COR2'
			'h1'		: tel = 'HI_1'
			'h2'		: tel = 'HI_2'
		ENDCASE
		XYOUTS, [.175], [.70 - tele / 20.], tel, SIZE = 1.5, COLOR = 0
	ENDFOR
;;	; Symbol in key for least squares line.
;;	OPLOT, [.150,.18], [.80,.80], LINESTYLE = 1, COLOR = 0
;;	XYOUTS, [.18], [.80], ' linear fit', SIZE = 1.5, COLOR = 0

	;-------------------------------------------------------
	; Calculate velocity of least squares line and print in 
	; key of Heigh-Time plot.
	; Velocity values >=100 resolve as integers.
	; Velocity values 10-99 resolve with 1 decimal place.
	; Velocity values 1-9 resolve with 2 decimal places.
	; Velocity values <1 resolve as floats.
	;-------------------------------------------------------
	elapsed_seconds = ( time[nfeatures-1] - time[0] ) * 24 * 3600
	slope = (REGRESS( (time - time[0]), height, CONST=offset, /DOUBLE ))[0]
	distance_km = (( time[nfeatures-1] * slope + offset ) - ( time[0] * slope + offset )) * km_sun_radius
	km_per_second = distance_km / elapsed_seconds
	CASE 1 OF
		( km_per_second GE 100 ) : km_per_second = STRTRIM( ROUND( km_per_second ), 2 )
		( km_per_second GE   1 ) : km_per_second = STRMID( STRTRIM( km_per_second, 2 ), 0, 4 )
		( km_per_second GE   1 ) : km_per_second = STRMID( STRTRIM( km_per_second, 2 ), 0, 4 )
		ELSE : km_per_second = STRTRIM( FLOAT( km_per_second ), 2 )
	ENDCASE
	XYOUTS, [.150], [.90], km_per_second + ' km/second', SIZE = 1.5, COLOR = 0
	WSET,curwin
ENDIF
IF (import_raw_xyz EQ 0) THEN BEGIN
	;-------------------------------------------------------
	; Restrict lat and lon to 2 decimal places.
	;-------------------------------------------------------
	CASE 1 OF
		( lat_mean GE 100 ) : lat_mean = STRMID(STRTRIM(lat_mean,2),0,6)
		( lat_mean GE 10 ) : lat_mean = STRMID(STRTRIM(lat_mean,2),0,5)
		( lat_mean GE 0 ) : lat_mean = STRMID(STRTRIM(lat_mean,2),0,4)
		( lat_mean GT -10 ) : lat_mean = STRMID(STRTRIM(lat_mean,2),0,5)
		( lat_mean GT -100 ) : lat_mean = STRMID(STRTRIM(lat_mean,2),0,6)
		ELSE : lat_mean = STRMID(STRTRIM(lat_mean,2),0,7)
	ENDCASE
	CASE 1 OF
		( lon_mean GE 100 ) : lon_mean = STRMID(STRTRIM(lon_mean,2),0,6)
		( lon_mean GE 10 ) : lon_mean = STRMID(STRTRIM(lon_mean,2),0,5)
		( lon_mean GE 0 ) : lon_mean = STRMID(STRTRIM(lon_mean,2),0,4)
		( lon_mean GT -10 ) : lon_mean = STRMID(STRTRIM(lon_mean,2),0,5)
		( lon_mean GT -100 ) : lon_mean = STRMID(STRTRIM(lon_mean,2),0,6)
		ELSE : lon_mean = STRMID(STRTRIM(lon_mean,2),0,7)
	ENDCASE
	CASE 1 OF
	        ( lat_mean_rmax GE 100 ) : lat_mean_rmax = STRMID(STRTRIM(lat_mean_rmax,2),0,6)
	        ( lat_mean_rmax GE 10 ) : lat_mean_rmax = STRMID(STRTRIM(lat_mean_rmax,2),0,5)
	        ( lat_mean_rmax GE 0 ) : lat_mean_rmax = STRMID(STRTRIM(lat_mean_rmax,2),0,4)
	        ( lat_mean_rmax GT -10 ) : lat_mean_rmax = STRMID(STRTRIM(lat_mean_rmax,2),0,5)
	        ( lat_mean_rmax GT -100 ) : lat_mean_rmax = STRMID(STRTRIM(lat_mean_rmax,2),0,6)
	        ELSE : lat_mean_rmax = STRMID(STRTRIM(lat_mean_rmax,2),0,7)
	ENDCASE
	CASE 1 OF
	        ( lon_mean_rmax GE 100 ) : lon_mean_rmax = STRMID(STRTRIM(lon_mean_rmax,2),0,6)
	        ( lon_mean_rmax GE 10 ) : lon_mean_rmax = STRMID(STRTRIM(lon_mean_rmax,2),0,5)
	        ( lon_mean_rmax GE 0 ) : lon_mean_rmax = STRMID(STRTRIM(lon_mean_rmax,2),0,4)
	        ( lon_mean_rmax GT -10 ) : lon_mean_rmax = STRMID(STRTRIM(lon_mean_rmax,2),0,5)
	        ( lon_mean_rmax GT -100 ) : lon_mean_rmax = STRMID(STRTRIM(lon_mean_rmax,2),0,6)
	        ELSE : lon_mean_rmax = STRMID(STRTRIM(lon_mean_rmax,2),0,7)
	ENDCASE
	CASE 1 OF
	        ( lon_sd_rmax GE 100 ) : lon_sd_rmax = STRMID(STRTRIM(lon_sd_rmax,2),0,6)
	        ( lon_sd_rmax GE 10 ) : lon_sd_rmax = STRMID(STRTRIM(lon_sd_rmax,2),0,5)
	        ( lon_sd_rmax GE 0 ) : lon_sd_rmax = STRMID(STRTRIM(lon_sd_rmax,2),0,4)
	        ( lon_sd_rmax GT -10 ) : lon_sd_rmax = STRMID(STRTRIM(lon_sd_rmax,2),0,5)
	        ( lon_sd_rmax GT -100 ) : lon_sd_rmax = STRMID(STRTRIM(lon_sd_rmax,2),0,6)
	        ELSE : lon_sd_rmax = STRMID(STRTRIM(lon_sd_rmax,2),0,7)
	ENDCASE
	CASE 1 OF
	        ( lat_sd_rmax GE 100 ) : lat_sd_rmax = STRMID(STRTRIM(lat_sd_rmax,2),0,6)
	        ( lat_sd_rmax GE 10 ) : lat_sd_rmax = STRMID(STRTRIM(lat_sd_rmax,2),0,5)
	        ( lat_sd_rmax GE 0 ) : lat_sd_rmax = STRMID(STRTRIM(lat_sd_rmax,2),0,4)
	        ( lat_sd_rmax GT -10 ) : lat_sd_rmax = STRMID(STRTRIM(lat_sd_rmax,2),0,5)
	        ( lat_sd_rmax GT -100 ) : lat_sd_rmax = STRMID(STRTRIM(lat_sd_rmax,2),0,6)
	        ELSE : lat_sd_rmax = STRMID(STRTRIM(lat_sd_rmax,2),0,7)
	ENDCASE
	;lat_mean = (STRSPLIT(STRTRIM(lat_mean,2),'.',/EXTRACT))[0] + '.' + STRMID((STRSPLIT(lat_mean,'.',/EXTRACT))[1],0,2)
	;lon_mean = (STRSPLIT(STRTRIM(lon_mean,2),'.',/EXTRACT))[0] + '.' + STRMID((STRSPLIT(lon_mean,'.',/EXTRACT))[1],0,2)
	IF (nfeatures GT 1) THEN BEGIN
		XYOUTS,[.150],[.80],system[0]+' Lon = ' + STRTRIM(lon_mean,2) + ' (' + STRTRIM(lon_mean_rmax,2) + ' +/- '+STRTRIM(lon_sd_rmax,2)+')',SIZE=1.5,COLOR=0
		XYOUTS,[.150],[.75],system[0]+' Lat = ' + STRTRIM(lat_mean,2) + ' (' + STRTRIM(lat_mean_rmax,2) + ' +/- '+STRTRIM(lat_sd_rmax,2)+')',SIZE=1.5,COLOR=0

		current_datetime_string = get_current_datetime_string()
		scc_ht_plot_filename = 'scc_ht_plot' + '_saved_' + current_datetime_string + '.png'
		WRITE_PNG,scc_ht_plot_filename,TVRD()

		;-------------------------------------------------------
		; Write txt file:
		; #Velocity = your linear fit
		; #t_lin_rs = t0, linear extrapolation of time to center of sun
		; #[Carrington|Stonyhurst] Lat (Lat max), lat mean (lat of highest feature)
		; #[Carrington|Stonyhurst] Lon (Lon max), lon mean (lon of highest feature)
		; # HEIGHT   DATE     TIME     TEL  xyzFilename
		;   2.93 2007/05/19 13:24:04    EUVI  filename.xyz
		;   4.33 2007/05/19 13:36:04   COR1  filename.xyz
		;-------------------------------------------------------
		scc_ht_info_filename = 'scc_ht_info' + '_saved_' + current_datetime_string + '.txt'
		OPENW, lun, scc_ht_info_filename, /GET_LUN, WIDTH = 160
		IF (nfeatures GT 1) THEN BEGIN
			PRINTF, lun, '#Velocity = ', km_per_second + ' km/second'
			PRINTF, lun, '#t_lin_rs = ', tzero_datetime
			PRINTF, lun, '#'+system[0]+' Lon = ' + STRTRIM(lon_mean,2) + ' (' + STRTRIM(lon_mean_rmax,2) + ' +/- '+STRTRIM(lon_sd_rmax,2)+')'
			PRINTF, lun, '#'+system[0]+' Lat = ' + STRTRIM(lat_mean,2) + ' (' + STRTRIM(lat_mean_rmax,2) + ' +/- '+STRTRIM(lat_sd_rmax,2)+')'
		ENDIF
		PRINTF, lun, '#'  + STRING(9B) + 'HEIGHT' + STRING(9B) + STRING(9B) + 'DATE' + STRING(9B) + STRING(9B) + 'TIME' + STRING(9B) + STRING(9B) + 'TEL' + STRING(9B) + STRING(9B) + STRING(9B) + 'xyzFilename'
		FOR i = 0, N_ELEMENTS(height) - 1 DO BEGIN
			hght = STRTRIM( height[i], 2 )
			datestring =	STRTRIM( FIX( STRMID( date[i], 4, 2 ) ), 2 )+ $
							'/'+STRTRIM( FIX( STRMID( date[i], 6, 2 ) ), 2 ) + $
							'/' + STRMID( date[i], 0, 4 )
			timestring =	STRMID( date[i], 9, 2 ) + $
							':' + STRMID( date[i], 11, 2 ) + $
							':' + STRMID( date[i], 13, 2 )
			CASE telescope[i] OF
				'171'	: tel = 'EUVI 171'
				'195'	: tel = 'EUVI 195'
				'284'	: tel = 'EUVI 284'
				'304'	: tel = 'EUVI 304'
				'eu 171'	: tel = 'EUVI 171'
				'eu 195'	: tel = 'EUVI 195'
				'eu 284'	: tel = 'EUVI 284'
				'eu 304'	: tel = 'EUVI 304'
				'c1'		: tel = 'COR1'
				'c2'		: tel = 'COR2'
				'h1'		: tel = 'HI_1'
				'h2'		: tel = 'HI_2'
			ENDCASE
			this_xyz_filename = FILE_BASENAME( xyz_filenames[i] )
			IF (STRMID(tel,0,4) EQ 'EUVI') THEN BEGIN
				PRINTF, lun, STRING(9B), hght, STRING(9B), datestring, STRING(9B), timestring, STRING(9B), tel, STRING(9B), this_xyz_filename
			ENDIF ELSE BEGIN
				PRINTF, lun, STRING(9B), hght, STRING(9B), datestring, STRING(9B), timestring, STRING(9B), tel, STRING(9B), STRING(9B), this_xyz_filename
			ENDELSE
		ENDFOR
		FREE_LUN,lun

	ENDIF
ENDIF

;-------------------------------------------------------
; Call LINE3D viewing widget.
;-------------------------------------------------------
IF (import_raw_xyz EQ 0) THEN BEGIN
	;-------------------------------------------------------
	; Normalize xyz to unit sphere.
	;-------------------------------------------------------
	xyz = xyz / solar_radius_km
	LINE3D, xyz, xyz_filenames[nfeatures-1], tlb, TIMESTEPS = date, XYZ_MEAN = [x_mean,y_mean,z_mean]
ENDIF ELSE BEGIN
	LINE3D, xyz, tlb, TIMESTEPS = date, SYSTEM = system[0], XYZ_MEAN = [x_mean,y_mean,z_mean]
ENDELSE

IF (nfeatures GT 1) THEN BEGIN
	IF ( !D.WINDOW GE 0 ) THEN BEGIN
		WSHOW
	ENDIF
ENDIF
END