; This function was translated from Joan Burkepile's CME trajectory
; fitting program, fit.f.  This function takes an array of strucutres defined
; in df.pro and db.pro.  The structure used in the array of structures is 
; defined in cmedata__define.pro.
;
; 2007-7-27 sitongia: correct calculation of polyfit weights and switch
;           to poly_fit

FUNCTION df_fit,data,fit

  LOADCT,13  ;Rainbow color table to use for plots
  TVLCT,r,g,b,/GET
  r[255]=255
  g[255]=255
  b[255]=255
  TVLCT,r,g,b

  ;--- Create data structure to return computed values in

  ret1={ret1,loop_start:'',cav_start:'',core_start:'',loop_v:0.0,cav_v:0.0,$
	core_v:0.0}
  ret2={ret2,loop_start:'',cav_start:'',core_start:'',loop_v:0.0,cav_v:0.0,$
	core_v:0.0,loop_h:0.0,cav_h:0.0,core_h:0.0,loop_accel:0.0,$
	cav_accel:0.0,core_accel:0.0}
  ;
  ; --- Determine what data is available
  ;
  loop_test=WHERE(data.feature EQ 'LOOP')
  cav_test=WHERE(data.feature EQ 'CAVITY')
  core_test=WHERE(data.feature EQ 'CORE')
  mk3_test=WHERE(data.dataset EQ 'MK3')
  c2_test=WHERE(data.dataset EQ 'C2')
  c3_test=WHERE(data.dataset EQ 'C3')
  IF (mk3_test[0] NE -1) THEN dset1=' MK3 ' ELSE dset1=''
  IF (c2_test[0] NE -1) THEN dset2=' C2 ' ELSE dset2=''
  IF (c3_test[0] NE -1) THEN dset3=' C3 ' ELSE dset3=''

  ;
  ; --- Open a .out file to put all of our information in
  ;
  dates=data.date_str[0]
  date=dates[0]
  doys=data.doy[0]
  doy=doys[0]
  tmp=STR_SEP(date,'/')
  yr=tmp[2]
  yrddoy=yr+'d'+STRCOMPRESS(STRING(doy,FORMAT="(i3.3)"),/REMOVE_ALL)
  ;path='/jira/d/tdarnell/outfiles/'
  SPAWN,'who am i',nametmp
  nm=str_sep(nametmp[0],' ')
  nm=strcompress(nm[0],/remove_all)
  filenme=nm+'.'+yrddoy+'.fit'+STRCOMPRESS(STRING(fit),/REMOVE_ALL)+'.out'
  ;OPENW,unit,path+filenme,/GET_LUN
  OPENW,unit,filenme,/GET_LUN
  PRINTF,unit,'This file generated by df_fit.pro'
  PRINTF,unit,'YRDDOY:  '+yrddoy
  PRINTF,unit,'Event date:  '+date
  PRINTF,unit,'Order Fit: '+STRCOMPRESS(STRING(fit))
  PRINTF,unit,''

  ;
  ; --- Perform fit on loop data
  ;
  IF (loop_test[0] NE -1) THEN BEGIN  ;{ 
    loop_data=data[WHERE(data.feature EQ 'LOOP')]
    tmp=SIZE(loop_data)
    endpoint=tmp[1]-1
    t_loop=loop_data[WHERE(loop_data.mins GE 0.0)].mins
    y_loop=loop_data[WHERE(loop_data.radius NE 0.0)].radius
    loop_err=loop_data[WHERE(loop_data.error NE 0.0)].error
    mins=t_loop
    loop_mins=mins
    t_loop[*]=t_loop[*]-t_loop[0]
    sz=SIZE(t_loop)
    loop_sz=sz
    ;x=FINDGEN(t_loop[sz[1]-1])
    ;x=FINDGEN(MAX(t_loop))
    ;loop_weights=SQRT(loop_data[WHERE(loop_data.error NE 0.0)].error)
    loop_errors=loop_data[WHERE(loop_data.error NE 0.0)].error

    ;
    ; Linear Fit =========================================================
    ;
    IF (fit EQ 1) THEN BEGIN  ;{
    ;loop_c=POLYFITW(t_loop,y_loop,loop_weights,1,yfit,yband,sigma,corr)
    loop_c=POLY_FIT(t_loop,y_loop,1,MEASURE_ERRORS=loop_errors,SIGMA=sigma)
    l1_start=mins[0]+(1.0-loop_c(0))/loop_c(1)     ;1st order start time (mins)
    ;l1_start=ABS(l1_start)
    lhr=LONG(l1_start/60)		    	   ;get the hours
    lmn=LONG(l1_start-60*lhr)			   ;get the minutes
    lsc=(l1_start - fix(l1_start)) * 60
    l1_str=STRING(lhr,lmn,lsc,FORMAT="(i2.2,':',i2.2,':',i2.2)")
    loop_vel=loop_c(1)*(6.96e5)/60		   ;Radial velocity
    lin_offset=ABS(mins[0]-l1_start)

    loop_subtitle2='!CLoop (Black):!C' + $
    'Radial Velocity: '+STRING(loop_vel,FORMAT="(F6.1)")+' km/s!C' + $
    'Start Time:      '+l1_str+' UT!C'
    PRINTF,unit,'Returned coefficients: '
    PRINTF,unit,'Constant term:  '+STRING(loop_c[0]*6.96e5)+' km OR ',loop_c[0]
    PRINTF,unit,'Linear term:    '+STRING(loop_c[1]*6.96e5)+' km OR ', $
		loop_c[1]
    PRINTF,unit,'Error in coeffs: '+STRING(sigma*6.96e5)+' km OR ',sigma
    PRINTF,UNIT,''
    PRINTF,unit,'Start Time:  '+l1_str
    PRINTF,unit,'Radial Velocity:  '+STRING(loop_vel)+' km/sec'
    PRINTF,unit,''
    ENDIF  ;}

    ;
    ;Quadratic fit ======================================================== 
    ;
    IF (fit EQ 2) THEN BEGIN  ;{
    ;loop_coeffs=POLYFITW(t_loop,y_loop,loop_weights,2,yfit,yband,sigma,corr)
    loop_coeffs=POLY_FIT(t_loop,y_loop,2,MEASURE_ERRORS=loop_errors,SIGMA=sigma)
    ; get start time and height
    c1=loop_coeffs(0)
    c2=loop_coeffs(1)
    c3=loop_coeffs(2)
    ; if quadratic solution is imaginary, then parabola does not intersect
    ; radius=1, set derivative = 1 and find minimum.
    IF (c2^2 LT (4*c3*(c1-1.0))) THEN BEGIN
      loop_tstart=-c2/(2*c3)
    ENDIF ELSE BEGIN
    ; if quadratic is real, then parabola intersects radius=1 and this is our
    ; start time for the event.
      quad1=-c2/(2*c3)
      quad2=SQRT(c2^2-4*c3*(c1-1.0))/(2*c3)
      loop_tstart=quad1+quad2
    ENDELSE
    l2_start=mins[0]+loop_tstart
    qhr=LONG(l2_start)/60
    qmn=LONG(l2_start-(60*qhr))
    qsc=(l2_start - fix(l2_start)) * 60
    l2_str=STRING(qhr,qmn,qsc,FORMAT="(i2.2,':',i2.2,':',i2.2)")

    ; start height is the quadratic evaluated at start time:
    loop_hstart=c1+(c2*loop_tstart)+(c3*loop_tstart^2)
    loop_accel=2*loop_coeffs(2)*(6.96e5)/3600
    ;get velocity at t final to get an idea of how fast it is going
    loop_finalv=(c2+c3*t_loop[sz[1]-1])*(6.96e5/60)
    quad_offset=ABS(loop_tstart)

    loop_subtitle='!CLoop (Black):!C'+ $
    'Radial Accel: '+STRING(loop_accel,FORMAT="(F6.3)")+' km/s!U2!N!C'+ $
    'Final speed: '+STRING(loop_finalv,FORMAT="(F6.1)")+' km/s!C' + $
    'Start Height: '+STRING(loop_hstart,FORMAT="(F3.1)")+' Solar Radii!C'+ $
    'Start Time: '+STRING(qhr,qmn,qsc,FORMAT="(i2.2,':',i2.2,':',i2.2)")+' UT!C'

    PRINTF,unit,'Returned Coefficients:'
    PRINTF,unit,'Constant term:  '+STRING(c1*6.96e5)+' km OR ',c1
    PRINTF,unit,'Linear term:    '+STRING(c2*6.96e5)+ 'km OR ',c2
    PRINTF,unit,'Quadratic term: '+STRING(c3*6.96e5)+' km OR ',c3
    PRINTF,unit,'Error in coeffs: '+STRING(sigma*6.96e5)+' km OR ',sigma
    ENDIF   ;}

    ;
    ; Cubic Fit =============================================================
    ;
    IF (fit EQ 3) THEN BEGIN  ;{
      ;3rd order fit
      ;third_coeffs=POLYFITW(t_loop,y_loop,loop_weights,3,yfit,yband,sigma,corr)
      third_coeffs=POLY_FIT(t_loop,y_loop,3,MEASURE_ERRORS=loop_errors,SIGMA=sigma)
      loop_coeffs=third_coeffs
      c1=third_coeffs(0)
      c2=third_coeffs(1)
      c3=third_coeffs(2)
      c4=third_coeffs(3)

    PRINTF,unit,'*********************** LOOP *****************************'
    PRINTF,unit,'Returned Coefficents:'
    PRINTF,unit,'Constant term:     '+STRING(c1*6.96e5)+' km OR ',c1
    PRINTF,unit,'Linear term:       '+STRING(c2*6.96e5)+' km OR ',c2
    PRINTF,unit,'Quadratic term:    '+STRING(c3*6.96e5)+' km OR ',c3
    PRINTF,unit,'Cubic term:        '+STRING(c4*6.96e5)+' km OR ',c4
    PRINTF,unit,'Error in coeffs:   '+STRING(sigma*6.96e5)+' km OR ',sigma
    PRINTF,unit,''

    ;start time is returned by df_starttimes.pro using the cubic 
    loop_tstart=df_starttimes(third_coeffs)
    tmp=WHERE(loop_tstart NE -1)
    tmp=SIZE(tmp)
    num_times=tmp[1]
    IF (num_times GT 1) THEN BEGIN
      l2_start=mins[0]+loop_tstart[0]
      l2_start2=mins[0]+loop_tstart[1] ; here is the second start time
      l2_str=format_time(l2_start)
      l2_str2=format_time(l2_start2)
      ; start height is the cubic evaluated at start time:
      loop_hstart=c1+(c2*loop_tstart[0])+(c3*loop_tstart[0]^2)+ $
		(c4*loop_tstart[0]^3)
      loop_hstart2=c1+(c2*loop_tstart[1])+(c3*loop_tstart[1]^2)+ $
		(c4*loop_tstart[1]^3)
      v1=c2+((2*c3)*loop_tstart[0])+(3*c4*loop_tstart[0]^2)
      v1=v1*6.96e5/60
      v2=c2+((2*c3)*loop_tstart[1])+(3*c4*loop_tstart[1]^2)
      v2=v2*6.96e5/60
      a1=((2*c3)+(6*c4*loop_tstart[0]))*6.96e5/3600
      a2=((2*c3)+(6*c4*loop_tstart[1]))*6.96e5/3600

      t=[l2_start,l2_start2]
      h=[loop_hstart,loop_hstart2]
      v=[v1,v2]
      a=[a1,a2]

      PRINTF,unit,'========== Start Time 1 ===========
      PRINTF,unit,''
      PRINTF,unit,'First start time:  '+l2_str
      PRINTF,unit,''
      PRINTF,unit,'Height:  ',loop_hstart
      PRINTF,unit,'Velocity: ',v1
      PRINTF,unit,'Accel:    ',a1
      PRINTF,unit,''
      PRINTF,unit,'========== Start time 2 ===========
      PRINTF,unit,''
      PRINTF,unit,'Second Start time:  '+l2_str2
      PRINTF,unit,''
      PRINTF,unit,'Height:  ',loop_hstart2
      PRINTF,unit,'Velocity: ',v2
      PRINTF,unit,'Accel:    ',a2
      PRINTF,unit,''
    ENDIF ELSE BEGIN
      l2_start=mins[0]+loop_tstart[0]
      l2_str=format_time(l2_start)
       ; start height is the cubic evaluated at start time:
      loop_hstart=c1+(c2*loop_tstart[0])+(c3*loop_tstart[0]^2)+ $
		(c4*loop_tstart[0]^3)
      v1=c2+((2*c3)*loop_tstart[0])+(3*c4*loop_tstart[0]^2)
      v1=v1*6.96e5/60
      a1=((2*c3)+(6*c4*loop_tstart[0]))*6.96e5/3600

      t=[l2_start]
      h=[loop_hstart]
      v=[v1]
      a=[a1]

      PRINTF,unit,'========== Start Time 1 ===========
      PRINTF,unit,''
      PRINTF,unit,'First start time:  '+l2_str
      PRINTF,unit,''
      PRINTF,unit,'Height:  ',loop_hstart
      PRINTF,unit,'Velocity: ',v1
      PRINTF,unit,'Accel:    ',a1
      PRINTF,unit,''
    ENDELSE


    h1=c1+(c2*t_loop[0])+(c3*t_loop[0]^2)+(c4*t_loop[0]^3)
    h2=c1+(c2*t_loop[endpoint])+(c3*t_loop[endpoint]^2)+(c4*t_loop[endpoint]^3)

    ;accel1 is initial acceleration evaluated at first datapoint
    ;accel2 is the final acceleration evaluated at the last datapoint
    loop_accel1=((2*c3)+(6*c4*t_loop[0]))*6.96e5/3600
    loop_accel2=((2*c3)+(6*c4*t_loop[endpoint]))*6.96e5/3600

    ;get velocity at t0 to get an idea of how fast it is going
    ;t0 is the point where the acceleration is zero, where velocity is at a
    ;max or min, and the inflection point of the height/time plot.
    t0=-(2*c3)/(6*c4)
    h_tau=c1+(c2*t0)+(c3*t0^2)+(c4*t0^3) ; height at tau
    a_tau=((2*c3)+(6*c4*t0))*6.96e5/3600 ; just to make sure it is zero
    ta0=mins[0]+t0
    IF (ta0 GT 1440) THEN ta0=ta0-1440  ;allow fo going past midnight.
    t0_hr=LONG(ta0)/60
    t0_min=LONG(ta0-(60*t0_hr))
    t0_sec=LONG(ta0-fix(ta0))*60
    t0_str=STRING(t0_hr,t0_min,t0_sec,FORMAT="(i2.2,':',i2.2,':',i2.2)")
    a0=2*c3
    v0=c2
    loop_finalv=c2+((2*c3)*t_loop[endpoint])+(3*c4*t_loop[endpoint]^2)
    loop_finalv=loop_finalv*6.96e5/60
    loop_initialv=(c2+((2*c3)*t_loop[0])+(3*c4*t_loop[0]^2)*6.96e5)/60
    loop_t0v=(v0+(a0*t0)-(a0/(2*t0))*t0^2)*6.96e5/60
    quad_offset=ABS(loop_tstart)

    loop_accel=loop_accel1

    loop_subtitle='!CLoop (Black):!C'+ $
    'Initial Accel: '+STRING(loop_accel1,FORMAT="(F6.3)")+' km/s!U2!N!C'+ $
    'Final Accel: '+STRING(loop_accel2,FORMAT="(F6.3)")+' km/s!U2!N!C'+ $
    'Vel @ a=0: '+STRING(loop_t0v,FORMAT="(F6.1)")+' km/s!C' + $
    'Start Height: '+STRING(loop_hstart,FORMAT="(F3.1)")+' Solar Radii!C'+ $
    'Accel=0 at: '+t0_str+' UT!C' + $
    'Start Time: '+l2_str

    PRINTF,unit,'======= AT FIRST DATAPOINT =========='
    PRINTF,unit,''
    PRINTF,unit,'Time:  ',format_time(mins[0])
    PRINTF,unit,'Height: ',h1
    PRINTF,unit,'Velocity:  ',loop_initialv
    PRINTF,unit,'accel:  ',loop_accel1
    PRINTF,unit,''
    PRINTF,unit,'======= AT LAST DATAPOINT =========='
    PRINTF,unit,''
    PRINTF,unit,'Time:  ',format_time(mins[endpoint])
    PRINTF,unit,'Height: ',h2
    PRINTF,unit,'Velocity:  ',loop_finalv
    PRINTF,unit,'accel:  ',loop_accel2
    PRINTF,unit,''
    PRINTF,unit,'======= AT TAU (a=0) =========='
    PRINTF,unit,''
    PRINTF,unit,'Time:  ',t0_str
    PRINTF,unit,'Height:  ',h_tau
    PRINTF,unit,'Velocity:  ',loop_t0v
    PRINTF,unit,'accel:  ',a_tau
    ENDIF  ;}

  ENDIF ELSE BEGIN
    loop_subtitle=''
    loop_subtitle2=''
    l1_str=''
    loop_vel=0.0
    loop_accel=0.0
    loop_accel1=0.0
    loop_accel2=0.0
    loop_finalv=0.0
    loop_hstart=0.0
    l2_str=''
  ENDELSE  ;}

  ;
  ; --- Perform fit on cavity data
  ;
  IF (cav_test[0] NE -1) THEN BEGIN
    cav_data=data[WHERE(data.feature EQ 'CAVITY')]
    t_cavity=cav_data[WHERE(cav_data.mins GE 0.0)].mins
    y_cavity=cav_data[WHERE(cav_data.radius NE 0.0)].radius
    tmp=SIZE(cav_data)
    endpoint=tmp[1]-1
    cav_mins=t_cavity
    t_cavity[*]=t_cavity[*]-t_cavity[0]
    cav_sz=SIZE(t_cavity)
    cav_x=FINDGEN(t_cavity[cav_sz[1]-1])
    cav_err=cav_data[WHERE(cav_data.error NE 0.0)].error
    ;cav_weights=SQRT(cav_data[WHERE(cav_data.error NE 0.0)].error)
    cav_errors=cav_data[WHERE(cav_data.error NE 0.0)].error

    ;
    ; Linear fit ==========================================================
    ;
    IF (fit EQ 1) THEN BEGIN  ;{
    ;1st order fit
    ;cav_c=POLYFITW(t_cavity,y_cavity,cav_weights,1)
    cav_c=POLY_FIT(t_cavity,y_cavity,1,MEASURE_ERRORS=cav_errors)
    c1_start=cav_mins[0]+(1.0-cav_c(0))/cav_c(1)   ;1st order start time (mins)
    ;c1_start=ABS(c1_start)
    chr=LONG(c1_start/60)		    	   ;get the hours
    cmn=LONG(c1_start-60*chr)			   ;get the minutes
    csc=(c1_start - fix(c1_start)) * 60
    c1_str=STRING(chr,cmn,csc,FORMAT="(i2.2,':',i2.2,':',i2.2)")
    cav_vel=cav_c(1)*(6.96e5)/60		   ;Radial velocity
    cav_subtitle2='!CCavity (Blue):!C' + $
      'Radial Velocity: '+STRING(cav_vel,FORMAT="(F5.1)")+' km/s!C' + $
      'Start Time:      '+c1_str+' UT!C'
    ENDIF  ;}

    ;
    ; Quadratic fit =======================================================
    ;
    IF (fit EQ 2) THEN BEGIN  ;{
    ;cav_coeffs=POLYFITW(t_cavity,y_cavity,cav_weights,2)
    cav_coeffs=POLY_FIT(t_cavity,y_cavity,2,MEASURE_ERRORS=cav_errors)
    ; get start time and height
    c1=cav_coeffs(0)
    c2=cav_coeffs(1)
    c3=cav_coeffs(2)
    IF (c2^2 LT (4*c3*(c1-1.0))) THEN BEGIN
      cav_tstart=-c2/(2*c3)
    ENDIF ELSE BEGIN
    ; if quadratic is real, then parabola intersects radius=1 and this is our
    ; start time for the event.
      quad1=-c2/(2*c3)
      quad2=SQRT(c2^2-4*c3*(c1-1.0))/(2*c3)
      cav_tstart=quad1+quad2
    ENDELSE
    c_start=cav_mins[0]+cav_tstart
    cav_str=format_time(c_start)
    c_hr=LONG(c_start/60)
    c_min=LONG(c_start-60*c_hr)
    c_sec=(c_start - fix(c_start)) * 60
    cav2_str=STRING(c_hr,c_min,c_sec,FORMAT="(i2.2,':',i2.2,':',i2.2)")

    ; start height is the quadratic evaluated at start time:
    cav_hstart=c1+c2*cav_tstart+c3*cav_tstart^2
    cav_accel=2*cav_coeffs(2)*(6.96e5)/3600
    ;get velocity at t-final to get an idea of how fast it is going
    cav_finalv=(c2+c3*t_cavity[cav_sz[1]-1])*(6.96e5/60)
    cavq_offset=ABS(cav_tstart)

    ;cav_vel=POLYFITW(t_cavity,y_cavity,cav_weights,1)
    cav_vel=POLY_FIT(t_cavity,y_cavity,1,MEASURE_ERRORS=cav_errors)
    cav_vel=cav_vel(1)*(6.96e5)/60

    cav_subtitle='!CCavity (Blue):!C'+ $
    'Radial Accel: '+STRING(cav_accel,FORMAT="(F6.3)")+' km/s!U2!N!C'+ $
    'Final Speed:  '+STRING(cav_finalv,FORMAT="(F6.1)")+' km/s!C'+ $
    'Start Height: '+STRING(cav_hstart,FORMAT="(F3.1)")+' Solar Radii!C'+ $
    'Start Time:   '+STRING(c_hr,c_min,c_sec, $
		FORMAT="(i2.2,':',i2.2,':',i2.2)")+' UT!C'
    ENDIF  ;}
    ;
    ; Cubic Fit =============================================================
    ;
    IF (fit EQ 3) THEN BEGIN  ;{
      ;3rd order fit
      ;third_coeffs=POLYFITW(t_cavity,y_cavity,cav_weights,3,yfit,yband,sigma,corr)
      third_coeffs=POLY_FIT(t_cavity,y_cavity,3,MEASURE_ERRORS=cav_errors,SIGMA=sigma)
      cav_coeffs=third_coeffs
      c1=third_coeffs(0)
      c2=third_coeffs(1)
      c3=third_coeffs(2)
      c4=third_coeffs(3)

    PRINTF,unit,'********************** CAVITY ****************************'
    PRINTF,unit,''
    PRINTF,unit,'Constant term:     '+STRING(c1*6.96e5)+' km OR ',c1
    PRINTF,unit,'Linear term:       '+STRING(c2*6.96e5)+' km OR ',c2
    PRINTF,unit,'Quadratic term:    '+STRING(c3*6.96e5)+' km OR ',c3
    PRINTF,unit,'Cubic term:        '+STRING(c4*6.96e5)+' km OR ',c4
    PRINTF,unit,'Error in coeffs:   '+STRING(sigma*6.96e5)+' km OR ',sigma
    PRINTF,unit,''

    ;start time is returned by df_starttimes.pro using the cubic
    cav_tstart=df_starttimes(cav_coeffs)
    tmp=WHERE(cav_tstart NE -1)
    tmp=SIZE(tmp)
    num_times=tmp[1]
    IF (num_times GT 1) THEN BEGIN
      c_start=cav_mins[0]+cav_tstart[0]
      c_start2=cav_mins[0]+cav_tstart[1]
      cav_str=format_time(c_start)
      cav_str2=format_time(c_start2)
      ; start height is the cubic evaluated at start time:
      cav_hstart =c1+(c2*cav_tstart[0])+(c3*cav_tstart[0]^2)+ $
		(c4*cav_tstart[0]^3)
      cav_hstart2 =c1+(c2*cav_tstart[1])+(c3*cav_tstart[1]^2)+ $
		(c4*cav_tstart[1]^3)
      v1=c2+((2*c3)*cav_tstart[0])+(3*c4*cav_tstart[0]^2)
      v1=v1*6.96e5/60
      v2=c2+((2*c3)*cav_tstart[1])+(3*c4*cav_tstart[1]^2)
      v2=v2*6.96e5/60
      a1=((2*c3)+(6*c4*cav_tstart[0]))*6.96e5/3600
      a2=((2*c3)+(6*c4*cav_tstart[1]))*6.96e5/3600

      tc=[c_start,c_start2]
      hc=[cav_hstart,cav_hstart2]
      vc=[v1,v2]
      ac=[a1,a2]

      PRINTF,unit,'========== Start Time 1 ===========
      PRINTF,unit,''
      PRINTF,unit,'First start time:  '+cav_str
      PRINTF,unit,''
      PRINTF,unit,'Height:  ',cav_hstart
      PRINTF,unit,'Velocity: ',v1
      PRINTF,unit,'Accel:    ',a1
      PRINTF,unit,''
      PRINTF,unit,'========== Start time 2 ===========
      PRINTF,unit,''
      PRINTF,unit,'Second Start time:  '+cav_str2
      PRINTF,unit,''
      PRINTF,unit,'Height:  ',cav_hstart2
      PRINTF,unit,'Velocity: ',v2
      PRINTF,unit,'Accel:    ',a2
      PRINTF,unit,''
    ENDIF ELSE BEGIN
      c_start=cav_mins[0]+cav_tstart[0]
      cav_str=format_time(c_start)
      ; start height is the cubic evaluated at start time:
      cav_hstart=c1+(c2*cav_tstart[0])+(c3*cav_tstart[0]^2)+ $
		(c4*cav_tstart[0]^3)
      v1=c2+((2*c3)*cav_tstart[0])+(3*c4*cav_tstart[0]^2)
      v1=v1*6.96e5/60
      a1=((2*c3)+(6*c4*cav_tstart[0]))*6.96e5/3600

      tc=[c_start]
      hc=[cav_hstart]
      vc=[v1]
      ac=[a1]

      PRINTF,unit,'========== Start Time 1 ===========
      PRINTF,unit,''
      PRINTF,unit,'First start time:  '+cav_str
      PRINTF,unit,''
      PRINTF,unit,'Height:  ',cav_hstart
      PRINTF,unit,'Velocity: ',v1
      PRINTF,unit,'Accel:    ',a1
      PRINTF,unit,''
    ENDELSE

    hc1=c1+(c2*t_cavity[0])+(c3*t_cavity[0]^2)+(c4*t_cavity[0]^3)
    hc2=c1+(c2*t_cavity[endpoint])+(c3*t_cavity[endpoint]^2)+ $
	(c4*t_cavity[endpoint]^3)

    ;accel1 is the initial acceleration
    ;accel2 is the final acceleration evaluated at the last datapoint
    cav_accel1=2*c3
    cav_accel2=((2*c3)+(6*c4*t_cavity[endpoint]))*6.96e5/3600

    ;get velocity at t0 to get an idea of how fast it is going
    ;t0 is the point where the acceleration is zero, where velocity is at a
    ;maximum, and the inflection point of the height/time plot.
    t0=-(2*c3)/(6*c4)
    h_tau=c1+(c2*t0)+(c3*t0^2)+(c4*t0^3)  ; height at tau
    a_tau=((2*c3)+(6*c4*t0))*6.96e5/3600  ; just to make sure it is zero
    ta0=cav_mins[0]+t0
    IF (ta0 GT 1440) THEN ta0=ta0-1440
    t0_hr=LONG(ta0)/60
    t0_min=LONG(ta0-(60*t0_hr))
    t0_sec=LONG(ta0-fix(ta0))*60
    t0_str=STRING(t0_hr,t0_min,t0_sec,FORMAT="(i2.2,':',i2.2,':',i2.2)")
    a0=2*c3
    v0=c2
    cav_finalv=(c2+((2*c3)*t_cavity[endpoint])+(3*c4*t_cavity[endpoint]^2))
    cav_finalv=cav_finalv*6.96e5/60
    cav_initialv=c2+(2*c3*t_cavity[0])+(3*c4*t_cavity[0]^2)
    cav_initialv=cav_initialv*6.96e5/60
    cav_t0v=(v0+(a0*t0)-(a0/(2*t0))*t0^2)*6.96e5/60
    cav_finalv=(c2+((2*c3)*t_cavity[cav_sz[1]-1])+ $
	(3*c4*t_cavity[cav_sz[1]-1]^2)*6.96e5)/60
    cavq_offset=ABS(cav_tstart)

    cav_accel=cav_accel1

    cav_subtitle='!CCavity (Blue):!C'+ $
    'Initial Accel: '+STRING(cav_accel1,FORMAT="(F6.3)")+' km/s!U2!N!C'+ $
    'Final Accel: '+STRING(cav_accel2,FORMAT="(F6.3)")+' km/s!U2!N!C'+ $
    'Max Velocity: '+STRING(cav_t0v,FORMAT="(F6.1)")+' km/s!C' + $
    'Start Height: '+STRING(cav_hstart,FORMAT="(F3.1)")+' Solar Radii!C'+ $
    'Accel = 0 at: '+t0_str+' UT!C' + $
    'Start Time: '+cav_str+' UT!C'

    PRINTF,unit,'======= AT FIRST DATAPOINT =========='
    PRINTF,unit,''
    PRINTF,unit,'Time:  ',format_time(cav_mins[0])
    PRINTF,unit,'Height: ',h1
    PRINTF,unit,'Velocity:  ',cav_initialv
    PRINTF,unit,'accel:  ',cav_accel1
    PRINTF,unit,''
    PRINTF,unit,'======= AT LAST DATAPOINT =========='
    PRINTF,unit,''
    PRINTF,unit,'Time:  ',format_time(cav_mins[endpoint])
    PRINTF,unit,'Height: ',h2
    PRINTF,unit,'Velocity:  ',cav_finalv
    PRINTF,unit,'accel:  ',cav_accel2
    PRINTF,unit,''
    PRINTF,unit,'======= AT TAU (a=0) =========='
    PRINTF,unit,''
    PRINTF,unit,'Time:  ',t0_str
    PRINTF,unit,'Height:  ',h_tau
    PRINTF,unit,'Velocity:  ',cav_t0v
    PRINTF,unit,'accel:  ',a_tau
    ENDIF  ;}

  ENDIF ELSE BEGIN
    cav_subtitle=''
    cav_subtitle2=''
    cav_vel=0.0
    cav_str=''
    c1_str=''
    cav_accel=0.0
    cav_finalv=0.0
    cav_hstart=0.0
    cav2_str=''
  ENDELSE

  ;
  ; --- Perform fit on core data
  ;
  IF (core_test[0] NE -1) THEN BEGIN
    core_data=data[WHERE(data.feature EQ 'CORE')]
    t_core=core_data[WHERE(core_data.mins GE 0.0)].mins
    y_core=core_data[where(core_data.radius NE 0.0)].radius
    tmp=SIZE(core_data)
    endpoint=tmp[1]-1
    core_mins=t_core
    t_core[*]=t_core[*]-t_core[0]
    core_sz=SIZE(t_core)
    core_x=FINDGEN(t_core[core_sz[1]-1])
    core_err=core_data[WHERE(core_data.error NE 0.0)].error
    ;core_weights=SQRT(core_data[WHERE(core_data.error NE 0.0)].error)
    core_errors=core_data[WHERE(core_data.error NE 0.0)].error

    ;
    ; Linear Fit ==========================================================
    ;
    IF (fit EQ 1) THEN BEGIN  ;{
    ;core_c=POLYFITW(t_core,y_core,core_weights,1)  ;1st order fit
    core_c=POLY_FIT(t_core,y_core,1,MEASURE_ERRORS=core_errors)
    c3_start=core_mins[0]+(1.0-core_c(0))/core_c(1) ;1st order start time (mins)
    ;c3_start=ABS(c3_start)
    whr=LONG(c3_start/60)		    	   ;get the hours
    wmn=LONG(c3_start-60*whr)			   ;get the minutes
    wsc=(c3_start - fix(c3_start)) * 60		   ;get the seconds
    c2_str=STRING(whr,wmn,wsc,FORMAT="(i2.2,':',i2.2,':',i2.2)")
    core_vel=core_c(1)*(6.96e5)/60		   ;Radial velocity
    core_subtitle2='!CCore (Red):!C' + $
      'Radial Velocity: '+STRING(core_vel,FORMAT="(F6.1)")+' km/s!C' + $
      'Start Time:      '+c2_str+' UT!C'
    ENDIF  ;}
    ;
    ; Quadratic Fit ========================================================
    ;
    IF (fit EQ 2) THEN BEGIN  ;{
    ;core_coeffs=POLYFITW(t_core,y_core,core_weights,2)
    core_coeffs=POLY_FIT(t_core,y_core,2,MEASURE_ERRORS=core_errors)
    ; get start time and height
    c1=core_coeffs(0)
    c2=core_coeffs(1)
    c3=core_coeffs(2)
    IF (c2^2 LT (4*c3*(c1-1.0))) THEN BEGIN
      core_tstart=-c2/(2*c3)
    ENDIF ELSE BEGIN
    ; if quadratic is real, then parabola intersects radius=1 and this is our
    ; start time for the event.
      quad1=-c2/(2*c3)
      quad2=SQRT(c2^2-4*c3*(c1-1.0))/(2*c3)
      core_tstart=quad1+quad2
    ENDELSE
    c2_start=core_mins[0]+core_tstart
    core_start_str=format_time(c2_start)
    c1_hr=LONG(c2_start/60)
    c1_min=LONG(c2_start-60*c1_hr)
    c1_sec=(c2_start - fix(c2_start)) * 60
    core_strt_str=STRING(c1_hr,c1_min,c1_sec,FORMAT="(i2.2,':',i2.2,':',i2.2)")

    ; start height is the quadratic evaluated at start time:
    core_hstart=c1+(c2*core_tstart)+(c3*core_tstart^2)

    ;core_vel=POLYFITW(t_core,y_core,core_weights,1)
    core_vel=POLY_FIT(t_core,y_core,1,MEASURE_ERRORS=core_errors)
    core_vel=core_vel(1)*(6.96e5)/60
    ;core_vel=core_coeffs(1)*(6.96e5)/60
    core_accel=2*core_coeffs(2)*(6.96e5)/3600
    core_finalv=(c2+c3*t_core[core_sz[1]-1])*(6.96e5/60)
    coreq_offset=ABS(core_tstart)

  core_subtitle='!CCore (Red):!C' + $
  'Radial Accel: '+STRING(core_accel,FORMAT="(F6.3)")+' km/s!U2!N!C'+ $
  'Final Speed:  '+STRING(core_finalv,FORMAT="(F6.1)")+' km/s!C'+ $
  'Start Height: '+STRING(core_hstart,FORMAT="(F3.1)")+' Solar Radii!C'+ $
  'Start Time:   '+STRING(c1_hr,c1_min,c1_sec, $
    FORMAT="(i2.2,':',i2.2,':',i2.2)") + ' UT'
  ENDIF  ;}
    ;
    ; Cubic Fit =============================================================
    ;
    IF (fit EQ 3) THEN BEGIN  ;{
      ;3rd order fit
      ;third_coeffs=POLYFITW(t_core,y_core,core_weights,3,yfit,yband,sigma,corr)
      third_coeffs=POLY_FIT(t_core,y_core,3,MEASURE_ERRORS=core_errors,SIGMA=sigma)
      core_coeffs=third_coeffs
      c1=third_coeffs(0)
      c2=third_coeffs(1)
      c3=third_coeffs(2)
      c4=third_coeffs(3)

    PRINTF,unit,'********************* CORE ****************************'
    PRINTF,unit,''
    PRINTF,unit,'Constant term:     '+STRING(c1*6.96e5)+' km OR ',c1
    PRINTF,unit,'Linear term:       '+STRING(c2*6.96e5)+' km OR ',c2
    PRINTF,unit,'Quadratic term:    '+STRING(c3*6.96e5)+' km OR ',c3
    PRINTF,unit,'Cubic term:        '+STRING(c4*6.96e5)+' km OR ',c4
    PRINTF,unit,'Error in coeffs:   '+STRING(sigma*6.96e5)+' km OR ',sigma
    PRINTF,unit,''

    ;start time is returned by df_starttimes.pro using the cubic
    core_tstart=df_starttimes(core_coeffs)
    tmp=WHERE(core_tstart NE -1)
    tmp=SIZE(tmp)
    num_times=tmp[1]		;get the number of start times returned.
    IF (num_times GT 1) THEN BEGIN
      c2_start=core_mins[0]+core_tstart[0]
      c2_start2=core_mins[0]+core_tstart[1]
      core_start_str=format_time(c2_start)
      core_start_str2=format_time(c2_start2)

      ; start height is the cubic evaluated at start time:
      core_hstart=c1+(c2*core_tstart[0])+(c3*core_tstart[0]^2)+ $
		(c4*core_tstart[0]^3)
      core_hstart2=c1+(c2*core_tstart[1])+(c3*core_tstart[1]^2)+ $
		(c4*core_tstart[1]^3)
      v1=c2+((2*c3)*core_tstart[0])+(3*c4*core_tstart[0]^2)
      v1=v1*6.96e5/60
      v2=c2+((2*c3)*core_tstart[1])+(3*c4*core_tstart[1]^2)
      v2=v2*6.96e5/60
      a1=((2*c3)+(6*c4*core_tstart[0]))*6.96e5/3600
      a2=((2*c3)+(6*c4*core_tstart[1]))*6.96e5/3600

      tcc=[c2_start,c2_start2]
      hcc=[core_hstart,core_hstart2]
      vcc=[v1,v2]
      acc=[a1,a2]

      PRINTF,unit,'========== Start Time 1 ===========
      PRINTF,unit,''
      PRINTF,unit,'First start time:  '+core_start_str
      PRINTF,unit,''
      PRINTF,unit,'Height:  ',core_hstart
      PRINTF,unit,'Velocity: ',v1
      PRINTF,unit,'Accel:    ',a1
      PRINTF,unit,''
      PRINTF,unit,'========== Start time 2 ===========
      PRINTF,unit,''
      PRINTF,unit,'Second Start time:  '+core_start_str2
      PRINTF,unit,''
      PRINTF,unit,'Height:  ',core_hstart2
      PRINTF,unit,'Velocity: ',v2
      PRINTF,unit,'Accel:    ',a2
      PRINTF,unit,''
    ENDIF ELSE BEGIN
      c2_start=core_mins[0]+core_tstart[0]
      core_start_str=format_time(c2_start)
      ; start height is the cubic evaluated at start time:
      core_hstart=c1+(c2*core_tstart[0])+(c3*core_tstart[0]^2)+ $
		(c4*core_tstart[0]^3)
      v1=c2+((2*c3)*core_tstart[0])+(3*c4*core_tstart[0]^2)
      v1=v1*6.96e5/60
      a1=((2*c3)+(6*c4*core_tstart[0]))*6.96e5/3600

      tcc=[c2_start]
      hcc=[core_hstart]
      vcc=[v1]
      acc=[a1]
      PRINTF,unit,'========== Start Time 1 ===========
      PRINTF,unit,''
      PRINTF,unit,'First start time:  '+core_start_str
      PRINTF,unit,''
      PRINTF,unit,'Height:  ',core_hstart
      PRINTF,unit,'Velocity: ',v1
      PRINTF,unit,'Accel:    ',a1
      PRINTF,unit,''
    ENDELSE

    hcc1=c1+(c2*t_core[0])+(c3*t_core[0]^2)+(c4*t_core[0]^3)
    hcc2=c1+(c2*t_core[endpoint])+(c3*t_core[endpoint]^2)+ $
	(c4*t_core[endpoint]^3)

    ;accel1 is initial acceleration 
    ;accel2 is the final acceleration evaluated at the last datapoint
    core_accel1=2*c3
    core_accel2=((2*c3)+(6*c4*t_core[endpoint]))*6.96e5/3600

    ;get velocity at t final to get an idea of how fast it is going
    ;t0 is the point where the acceleration is zero, where velocity is at a
    ;maximum, and the inflection point of the height/time plot.
    t0=-(2*c3)/(6*c4)
    h_tau=c1+(c2*t0)+(c3*t0^2)+(c4*t0^3)  ; height at tau
    a_tau=((2*c3)+(6*c4*t0))*6.96e5/3600  ; just to make sure it is zero
    ta0=core_mins[0]+t0
    IF (ta0 GT 1440) THEN ta0=ta0-1440
    t0_str=format_time(ta0)
    a0=2*c3
    v0=c2
    core_finalv=c2+(2*c3*t_core[endpoint])+(3*c4*t_core[endpoint]^2)
    core_finalv=core_finalv*6.96e5/60
    core_initialv=c2+(2*c3*t_core[0])+(3*c4*t_core[0]^2)
    core_initialv=core_finalv*6.96e5/60
    core_t0v=(v0+(a0*t0)-(a0/(2*t0))*t0^2)*6.96e5/60
    coreq_offset=ABS(core_tstart)

    core_accel=core_accel1

    core_subtitle='!CCore (Red):!C'+ $
    'Initial Accel: '+STRING(core_accel1,FORMAT="(F6.3)")+' km/s!U2!N!C'+ $
    'Final Accel: '+STRING(core_accel2,FORMAT="(F6.3)")+' km/s!U2!N!C'+ $
    'Max speed: '+STRING(core_t0v,FORMAT="(F6.1)")+' km/s!C' + $
    'Start Height: '+STRING(core_hstart,FORMAT="(F3.1)")+' Solar Radii!C'+ $
    'Accel = 0 at: '+t0_str+' UT!C' + $
    'Start Time: '+core_start_str+' UT!C'

    PRINTF,unit,'======= AT FIRST DATAPOINT =========='
    PRINTF,unit,''
    PRINTF,unit,'Time:  ',format_time(core_mins[0])
    PRINTF,unit,'Height: ',hcc1
    PRINTF,unit,'Velocity:  ',core_initialv
    PRINTF,unit,'accel:  ',core_accel1
    PRINTF,unit,''
    PRINTF,unit,'======= AT LAST DATAPOINT =========='
    PRINTF,unit,''
    PRINTF,unit,'Time:  ',format_time(core_mins[endpoint])
    PRINTF,unit,'Height: ',hcc2
    PRINTF,unit,'Velocity:  ',core_finalv
    PRINTF,unit,'accel:  ',core_accel2
    PRINTF,unit,''
    PRINTF,unit,'======= AT TAU (a=0) =========='
    PRINTF,unit,''
    PRINTF,unit,'Time:  ',t0_str
    PRINTF,unit,'Height:  ',h_tau
    PRINTF,unit,'Velocity:  ',core_t0v
    PRINTF,unit,'accel:  ',a_tau
    ENDIF  ;}

  ENDIF ELSE BEGIN
    core_subtitle=''
    core_subtitle2=''
    c2_str=''
    core_vel=0.0
    core_accel=0.0
    core_finalv=0.0
    core_hstart=0.0
    core_start_str=''
  ENDELSE

  !Y.RANGE(0)=1.0
  radii=data[WHERE(data.radius NE 0.0)].radius
  !Y.RANGE(1)=MAX(radii)

  mins=data[WHERE(data.mins GE 0.0)].mins
  sz=SIZE(mins)
  x=FINDGEN(mins[sz(1)-1]+5)
  x=x-mins[0]

  ;
  ; Linear Plot========================================================
  ;
  IF (fit EQ 1) THEN BEGIN  ;{
  IF (loop_test[0] NE -1) THEN orig=l1_start ELSE $
  IF (cav_test[0] NE -1) THEN orig=c1_start ELSE $
  orig=c3_start
  end_x=MAX(mins)

  !P.SUBTITLE=loop_subtitle2+cav_subtitle2+core_subtitle2

  PLOT,x,radii,BACKGROUND=255,COLOR=0, $
    TITLE='CME Trajectory Plot!C'+data[0].date_str+' - Linear Fit', $
    ;XRANGE=[mins[0]-ABS(orig),end_x], $
    XRANGE=[orig,end_x], $
    XTITLE='Time (UT)',YTITLE='Height (Solar Radii)',PSYM=1, $
    XTICKFORMAT='xticks'

  IF (loop_test[0] NE -1) THEN BEGIN
    t1_loop=t_loop+loop_mins[0]
    loop_x1=FINDGEN(t1_loop[loop_sz(1)-1]+5)
    loop_x1=loop_x1-loop_mins[0]
    func=loop_c(0)+(loop_c(1)*(loop_x1))
    loop_chi_squared=XSQ_TEST(y_loop,func[t1_loop])
    OPLOTERR,t1_loop,y_loop,loop_err,PSYM=1
    OPLOT,t1_loop,y_loop,COLOR=0,PSYM=1
    OPLOT,loop_c(0)+loop_c(1)*(loop_x1),COLOR=0
  ENDIF ELSE loop_chi_squared = -1

  IF (cav_test[0] NE -1) THEN BEGIN
    t1_cavity=t_cavity+cav_mins[0]
    cav_x1=FINDGEN(t1_cavity[cav_sz[1]-1]+5)
    cav_x1=cav_x1-cav_mins[0]
    func=cav_c(0)+(cav_c(1)*(cav_x1))
    cav_chi_squared=XSQ_TEST(y_cav,func[t1_cavity])
    OPLOTERR,t1_cavity,y_cavity,cav_err,PSYM=4
    OPLOT,t1_cavity,y_cavity,COLOR=60,PSYM=4
    OPLOT,cav_c(0)+cav_c(1)*cav_x1,COLOR=60
  ENDIF ELSE cav_chi_squared = -1

  IF (core_test[0] NE -1) THEN BEGIN
    t1_core=t_core+core_mins[0]
    core_x1=FINDGEN(t1_core[core_sz[1]-1]+5)
    core_x1=core_x1-core_mins[0]
    func=core_c(0)+(core_c(1)*(core_x1))
    core_chi_squared=XSQ_TEST(y_core,func[t1_core])
    OPLOTERR,t1_core,y_core,core_err,PSYM=4
    OPLOT,t1_core,y_core,COLOR=253,PSYM=4
    OPLOT,core_c(0)+core_c(1)*core_x1,COLOR=253
  ENDIF ELSE core_chi_squared = -1
  ENDIF  ;}


  ;
  ; Quadratic plot =====================================================
  ;
  IF (fit EQ 2) THEN BEGIN  ;{

  ; Get origin
    IF (loop_test[0] NE -1) THEN origin=l2_start ELSE $
    IF (cav_test[0] NE -1) THEN origin=c_start ELSE $
    origin=c2_start
  ;Get end of x-axis
    end_x=MAX(mins)
  ;
  ; Do annotation, (set !P.SUBTITLE)
  ;
   !P.SUBTITLE=loop_subtitle+cav_subtitle+core_subtitle 

  ;PLOT,t_loop,y_loop,BACKGROUND=255,COLOR=0, $
  PLOT,x,radii,BACKGROUND=255,COLOR=0, $
    TITLE='CME Trajectory Plot!C'+data[0].date_str+' - Quadratic Fit', $
    ;XRANGE=[mins[0]-ABS(origin),end_x], $
    XRANGE=[origin,end_x], $
    XTITLE='Time (UT)',YTITLE='Height (Solar Radii)',PSYM=1, $
    XTICKFORMAT='xticks'

  IF (loop_test[0] NE -1) THEN BEGIN
    t_loop=t_loop+loop_mins[0]
    loop_x=FINDGEN(t_loop[loop_sz(1)-1]+5)
    loop_x=loop_x-loop_mins[0]
    OPLOTERR,t_loop,y_loop,loop_err,PSYM=1
    OPLOT,t_loop,y_loop,COLOR=0,PSYM=1
    func=loop_coeffs(0)+(loop_coeffs(1)*(loop_x)) + $
    	(loop_coeffs(2)*(loop_x)^2)
    mn=WHERE(func EQ MIN(func))
    loop_chi_squared=XSQ_TEST(y_loop,func[t_loop])
    mask=func[0:mn[0]]  ; This will clear out unwanted part of parabola.
    OPLOT,func,COLOR=0
    OPLOT,mask,COLOR=255
    ;OPLOT,loop_coeffs(0)+(loop_coeffs(1)*(loop_x)) + $
    ;	(loop_coeffs(2)*(loop_x)^2),COLOR=0,MIN_VALUE=0.0
  ENDIF ELSE loop_chi_squared=-1

  IF (cav_test[0] NE -1) THEN BEGIN
    t_cavity=t_cavity+cav_mins[0]
    cav_x=FINDGEN(t_cavity[cav_sz[1]-1]+5)
    cav_x=cav_x-cav_mins[0]
    func=cav_coeffs(0)+cav_coeffs(1)*cav_x+cav_coeffs(2)*cav_x^2
    mn=WHERE(func EQ MIN(func))
    cav_chi_squared=XSQ_TEST(y_cavity,func[t_cavity])
    mask=func[0:mn[0]]  ; This will clear out unwanted part of parabola.
    OPLOTERR,t_cavity,y_cavity,cav_err,PSYM=4
    OPLOT,t_cavity,y_cavity,COLOR=60,PSYM=4
    OPLOT,func,COLOR=60
    OPLOT,mask,COLOR=255
    ;OPLOT,cav_coeffs(0)+cav_coeffs(1)*cav_x+cav_coeffs(2)*cav_x^2,COLOR=60
  ENDIF ELSE cav_chi_squared=-1

  IF (core_test[0] NE -1) THEN BEGIN
    t_core=t_core+core_mins[0]
    core_x=FINDGEN(t_core[core_sz[1]-1]+5)
    core_x=core_x-core_mins[0]
    func=core_coeffs(0)+core_coeffs(1)*core_x+core_coeffs(2)*core_x^2
    mn=WHERE(func EQ MIN(func))
    mask=func[0:mn[0]]  ; This will clear out unwanted part of parabola.
    core_chi_squared=XSQ_TEST(y_core,func[t_core])
    OPLOTERR,t_core,y_core,core_err,PSYM=2
    OPLOT,t_core,y_core,COLOR=253,PSYM=2
    OPLOT,func,COLOR=253
    OPLOT,mask,COLOR=255
    ;OPLOT,core_coeffs(0)+core_coeffs(1)*core_x+core_coeffs(2)*core_x^2, $
    ;  COLOR=253
  ENDIF ELSE core_chi_squared=-1
  ENDIF  ;}

  ;
  ; Cubic Plot ================================================================
  ;
  IF (fit EQ 3) THEN BEGIN  ;{

  ; Get origin
    IF (loop_test[0] NE -1) THEN origin=l2_start ELSE $
    IF (cav_test[0] NE -1) THEN origin=c_start ELSE $
    origin=c2_start
  ;Get end of x-axis
    end_x=MAX(mins)
  ;
  ; Do annotation, (set !P.SUBTITLE)
  ;
   !P.SUBTITLE=loop_subtitle+cav_subtitle+core_subtitle 

  ;
  ; Trajectory Plot
  ;
  ;PLOT,t_loop,y_loop,BACKGROUND=255,COLOR=0, $
  PLOT,x,radii,BACKGROUND=255,COLOR=0, $
    TITLE='CME Trajectory Plot!C'+data[0].date_str+' - Cubic Fit', $
    ;XRANGE=[mins[0]-ABS(origin),end_x], $
    XRANGE=[origin,end_x], $
    XTITLE='Time (UT)',YTITLE='Height (Solar Radii)',PSYM=1, $
    XTICKFORMAT='xticks'

  IF (loop_test[0] NE -1) THEN BEGIN
    t_loop=t_loop+loop_mins[0]
    loop_x=FINDGEN(t_loop[loop_sz(1)-1]+5)
    loop_x=loop_x-loop_mins[0]
    zero=loop_x
    zero[*]=0
    xaxis=FINDGEN(1000)+loop_mins[0]-400
    OPLOTERR,t_loop,y_loop,loop_err,PSYM=1
    OPLOT,t_loop,y_loop,COLOR=0,PSYM=1
    OPLOT,t,h,COLOR=250,PSYM=5  ;plot the three start times
    func=loop_coeffs(0)+(loop_coeffs(1)*(loop_x)) + $
    	(loop_coeffs(2)*(loop_x)^2+loop_coeffs(3)*(loop_x)^3)
    loop_vel_func=loop_coeffs(1)+(2*loop_coeffs(2)*loop_x)+ $
    		(3*loop_coeffs(3)*(loop_x)^2)
    loop_vel_func=loop_vel_func*6.96e5/60  ; put function in proper units
    loop_acc_func=2*loop_coeffs(2)+(6*loop_coeffs(3)*loop_x)
    loop_acc_func=loop_acc_func*6.96e5/3600
    mn=WHERE(func EQ MIN(func))
    loop_chi_squared=XSQ_TEST(y_loop,func[t_loop])
    ;print,'chi squared (loop): ',loop_chi_squared[0]
    OPLOT,func,COLOR=0
  ENDIF ELSE loop_chi_squared=-1

  IF (cav_test[0] NE -1) THEN BEGIN
    t_cavity=t_cavity+cav_mins[0]
    cav_x=FINDGEN(t_cavity[cav_sz[1]-1]+5)
    cav_x=cav_x-cav_mins[0]
    func=cav_coeffs(0)+cav_coeffs(1)*cav_x+cav_coeffs(2)*cav_x^2+ $
      cav_coeffs(3)*cav_x^3
    cav_vel_func=cav_coeffs(1)+(2*cav_coeffs(2)*cav_x)+ $
	(3*cav_coeffs(3)*(cav_x)^2)
    cav_vel_func=cav_vel_func*6.96e5/60
    cav_acc_func=2*cav_coeffs(2)+(6*cav_coeffs(3)*cav_x)
    cav_acc_func=cav_acc_func*6.96e5/3600
    mn=WHERE(func EQ MIN(func))
    cav_chi_squared=XSQ_TEST(y_cavity,func[t_cavity])
    ;print,'chi squared (cavity): ',cav_chi_squared[0]
    OPLOTERR,t_cavity,y_cavity,cav_err,PSYM=4
    OPLOT,t_cavity,y_cavity,COLOR=60,PSYM=4
    OPLOT,func,COLOR=60
    OPLOT,tc,hc,COLOR=250,PSYM=5  ;plot the two start times.
  ENDIF ELSE cav_chi_squared=-1

  IF (core_test[0] NE -1) THEN BEGIN
    t_core=t_core+core_mins[0]
    core_x=FINDGEN(t_core[core_sz[1]-1]+5)
    core_x=core_x-core_mins[0]
    func=core_coeffs(0)+core_coeffs(1)*core_x+core_coeffs(2)*core_x^2 + $
      core_coeffs(3)*core_x^3
    core_vel_func=core_coeffs(1)+(2*core_coeffs(2)*core_x)+ $
	(3*core_coeffs(3)*(core_x)^2)
    core_vel_func=core_vel_func*6.96e5/60
    core_acc_func=2*core_coeffs(2)+(6*core_coeffs(3)*core_x)
    core_acc_func=core_acc_func*6.96e5/3600
    core_chi_squared=XSQ_TEST(y_core,func[t_core])
    ;print,'chi squared (core): ',core_chi_squared[0]
    mn=WHERE(func EQ MIN(func))
    OPLOTERR,t_core,y_core,core_err,PSYM=2
    OPLOT,t_core,y_core,COLOR=253,PSYM=2
    OPLOT,func,COLOR=253
    OPLOT,tcc,hcc,COLOR=250,PSYM=5  ;plot the two start times.
  ENDIF ELSE core_chi_squared=-1

  XYOUTS,3000,16000,dset1+dset2+dset3,/DEVICE
  XYOUTS,3000,15000,'Chi squared (loop): '+STRING(loop_chi_squared[0]),/DEVICE
  XYOUTS,3000,14500,'Chi squared (cav): '+STRING(cav_chi_squared[0]),/DEVICE
  XYOUTS,3000,14000,'Chi squared (core): '+STRING(core_chi_squared[0]),/DEVICE
  ;
  ; Velocity plot  (Cubic)
  ;
  !Y.RANGE(0)=0
  !Y.RANGE(1)=900

  PLOT,x,radii,BACKGROUND=255,COLOR=0, $
    TITLE='CME Velocity Plot!C'+data[0].date_str+' - Cubic Fit', $
    XRANGE=[origin,end_x], $
    XTITLE='Time (UT)',YTITLE='Velocity (km/sec)',PSYM=1, $
    XTICKFORMAT='xticks'
  
IF (loop_test[0] NE -1) AND (cav_test[0] NE -1) AND (core_test[0] NE -1) THEN $
  BEGIN
    OPLOT,loop_vel_func,COLOR=0
    OPLOT,cav_vel_func,COLOR=60
    OPLOT,core_vel_func,COLOR=253
    OPLOT,t,v,COLOR=250,PSYM=5  ;plot the three start times (loop)
    OPLOT,tc,vc,COLOR=250,PSYM=5  ;plot the three start times (cav)
    OPLOT,tcc,vcc,COLOR=250,PSYM=5  ;plot the three start times (core)
  ENDIF ELSE IF (loop_test[0] NE -1) AND (cav_test[0] NE -1) THEN $
  BEGIN
    OPLOT,loop_vel_func,COLOR=0
    OPLOT,cav_vel_func,COLOR=60
    OPLOT,t,v,COLOR=250,PSYM=5  ;plot the three start times (loop)
    OPLOT,tc,vc,COLOR=250,PSYM=5  ;plot the three start times (cav
  ENDIF ELSE BEGIN  ; only one dataset is here 
    IF (loop_test[0] NE -1) THEN BEGIN
      OPLOT,loop_vel_func,COLOR=0 
      OPLOT,t,v,COLOR=250,PSYM=5  ;plot the three start times (loop)
    ENDIF ELSE IF (cav_test[0] NE -1) THEN BEGIN
      OPLOT,cav_vel_func,COLOR=60
      OPLOT,tc,vc,COLOR=250,PSYM=5  ;plot the three start times (cav)
    ENDIF ELSE IF (core_test[0] NE -1) THEN BEGIN
      OPLOT,core_vel_func,COLOR=253
      OPLOT,tcc,vcc,COLOR=250,PSYM=5  ;plot the three start times (core)
    ENDIF
  ENDELSE
  ;
  ; Acceleration Plot
  ;
    !Y.RANGE(0)=-0.5
    !Y.RANGE(1)=0.5

    PLOT,x,radii,BACKGROUND=255,COLOR=0, $
        TITLE='CME Acceleration Plot!C'+data[0].date_str+' - Cubic Fit', $
      XRANGE=[origin,end_x], $
      XTITLE='Time (UT)',YTITLE='Acceleration (km/sec!U2!N)',PSYM=1, $
      XTICKFORMAT='xticks'
IF (loop_test[0] NE -1) AND (cav_test[0] NE -1) AND (core_test[0] NE -1) THEN $
  BEGIN
    OPLOT,xaxis,zero,COLOR=0  ;make a line at a=0
    OPLOT,loop_acc_func,COLOR=0
    OPLOT,cav_acc_func,COLOR=60
    OPLOT,core_acc_func,COLOR=253
    OPLOT,t,a,COLOR=250,PSYM=5  ;plot the three start times (loop)
    OPLOT,tc,ac,COLOR=250,PSYM=5  ;plot the three start times (cav)
    OPLOT,tcc,acc,COLOR=250,PSYM=5  ;plot the three start times (core)
  ENDIF ELSE IF (loop_test[0] NE -1) AND (cav_test[0] NE -1) THEN $
  BEGIN
    OPLOT,loop_acc_func,COLOR=0
    OPLOT,cav_acc_func,COLOR=60
    OPLOT,t,a,COLOR=250,PSYM=5  ;plot the three start times (loop)
    OPLOT,tc,ac,COLOR=250,PSYM=5  ;plot the three start times (cav)
  ENDIF ELSE BEGIN
    IF (loop_test[0] NE -1) THEN BEGIN
      OPLOT,loop_acc_func,COLOR=0 
      OPLOT,t,a,COLOR=250,PSYM=5  ;plot the three start times (loop)
    ENDIF ELSE IF (cav_test[0] NE -1) THEN BEGIN
      OPLOT,cav_acc_func,COLOR=60
      OPLOT,tc,ac,COLOR=250,PSYM=5  ;plot the three start times (cav)
    ENDIF ELSE IF (core_test[0] NE -1) THEN BEGIN
      OPLOT,core_acc_func,COLOR=253
      OPLOT,tcc,acc,COLOR=250,PSYM=5  ;plot the three start times (core)
    ENDIF
  ENDELSE
    
  ENDIF  ;}

  XYOUTS,3000,16000,dset1+dset2+dset3,/DEVICE
  XYOUTS,3000,15000,'Chi squared (loop): '+STRING(loop_chi_squared[0]),/DEVICE
  XYOUTS,3000,14500,'Chi squared (cav): '+STRING(cav_chi_squared[0]),/DEVICE
  XYOUTS,3000,14000,'Chi squared (core): '+STRING(core_chi_squared[0]),/DEVICE
  CLOSE,unit

  mins=0.0
  x=0.0
  origin=0.0
  end_x=0.0

  IF (fit EQ 1) THEN BEGIN
    ret1.loop_start=l1_str
    ret1.cav_start=c1_str
    ret1.core_start=c2_str

    ret1.loop_v=loop_vel
    ret1.cav_v=cav_vel
    ret1.core_v=core_vel

    RETURN,ret1
  ENDIF ELSE BEGIN
    ret2.loop_start=l2_str
    ret2.cav_start=cav_str
    ret2.core_start=core_start_str

    ret2.loop_v=loop_finalv
    ret2.cav_v=cav_finalv
    ret2.core_v=core_finalv

    ret2.loop_h=loop_hstart
    ret2.cav_h=cav_hstart
    ret2.core_h=core_hstart

    ret2.loop_accel=loop_accel
    ret2.cav_accel=cav_accel
    ret2.core_accel=core_accel

    RETURN,ret2
  ENDELSE

END