;
;  This function performs a 3rd order, least-squares fit of a given set
;	of data.  The computed coefficients are returned and a variety of
;	output is printed in a .out file.
;
FUNCTION df_fit3,data,origin=origin,end_x=end_x,noannotate=noannotate

  num_meas=SIZE(data.dataset)
  num_meas=num_meas[1]
  num_features=1
  num_datasets=1
  index=sort(data.dataset)
  srt=data[index].dataset
  FOR i=1,num_meas-1 DO BEGIN
      ;IF (data[i].dataset NE data[i-1].dataset) THEN $
      ;	num_datasets=num_datasets+1
      IF (srt[i] NE srt[i-1]) THEN num_datasets=num_datasets+1
      IF (data[i].feature NE data[i-1].feature) THEN $
	num_features=num_features+1
  ENDFOR

  print,'num_features:  ',num_features
  print,'num_datasets:  ',num_datasets

  date=str_sep(data[0].date_str,'/')
  yrddoy=date[2]+'d'+STRTRIM(STRING(data[0].doy,FORMAT="(i3.3)"),1)
  mins=data.mins
  radii=data.radius
  minutes=FINDGEN(MAX(mins))
  sz=SIZE(minutes)


  features=STRARR(num_features)
  datasets=STRARR(num_datasets)
  j=0
  k=0
  FOR i=1,num_meas-1 DO BEGIN
    IF (data[i].feature NE data[i-1].feature) THEN BEGIN
      features[j]=data[i-1].feature
      j=j+1
    ENDIF
  features[j]=data[i].feature
    IF (srt[i] NE srt[i-1]) THEN BEGIN
      datasets[k]=srt[i-1]
      k=k+1
    ENDIF
  datasets[k]=srt[i]
  ENDFOR
  print,features
  print,datasets

  ;
  ; --- Open a .out file to put all of our information in
  ;
  path='/bail/d/tdarnell/outfiles/'
  filename=yrddoy+'.fit3.out'
  OPENW,unit,path+filename,/GET_LUN
  PRINTF,unit,''
  PRINTF,unit,'Order Fit:  3'
  PRINTF,unit,''

  ;
  ; Cubic Fit =========================================================
  ;
  coeff_arr=FLTARR(num_features,4)
  starttimes=FLTARR(num_features,2)  ;each feature may have two start times
  timestrarr=STRARR(num_features,3)  ;two start times and t0
  startheights=FLTARR(num_features,2) ;one for each start time
  ;
  ; Velocites contain the velocities at the two starttimes, final, and t0
  ; velocity for each feature.
  velocities=FLTARR(num_features,5) 
  accels=FLTARR(num_features,4)
  functions=FLTARR(num_features,sz(1))
  vel_funcs=functions
  acc_funcs=functions

  FOR i=0, num_features-1 DO BEGIN  ;{

  xdata=data[where(data.feature EQ features[i])].mins
  ydata=data[where(data.feature EQ features[i])].radius
  errors=SQRT(data[where(data.feature EQ features[i])].error)
  offset=xdata[0]
  endpoint=N_ELEMENTS(xdata)-1

  xdata=xdata-offset

  coeffs=POLYFITW(xdata,ydata,errors,3,yfit,yband,sigma,corr)

  ;xdata=xdata+offset

  ;
  ; Calculate start times
  ;
  ;starttimes[i,*]=df_starttimes(coeffs) + offset
  starttimes[i,*]=df_starttimes(coeffs)
  valid=where(starttimes NE -1)
  index=where(starttimes EQ -1)
  ;starttimes[valid]=starttimes[valid]+offset

  ;
  ; Calculate t0, the inflection point of the fit.  This indicates
  ;  where acceleration=0, the slope of the plot is at a max.
  ;
  t0=-(2*coeffs[2])/(6*coeffs[3])
  ta0=t0+offset

  IF (starttimes[i,0] NE -1) THEN $
  timestrarr[i,0]=format_time(starttimes[i,0]+offset) ELSE $
  timestrarr[i,0]='Not applicable'
  IF (starttimes[i,1] NE -1) THEN $
  timestrarr[i,1]=format_time(starttimes[i,1]+offset) ELSE $
  timestrarr[i,1]='Not applicable'
  timestrarr[i,2]=format_time(ta0)

  ;
  ;Start height is the cubic evaluated at the start time.
  ;
  ;starttimes[valid]=starttimes[valid]-offset

  minutes=minutes-offset
  startheights[i,0:1]=coeffs(0)+(coeffs(1)*starttimes[i,*])+ $
	(coeffs(2)*starttimes[i,*]^2)+(coeffs(3)*starttimes[i,*]^3)
  IF (index[0] NE -1) THEN startheights[index] = 0.0
  functions[i,*]=coeffs(0)+(coeffs(1)*minutes)+coeffs(2)*minutes^2+ $
		coeffs(3)*minutes^3
  vel_funcs[i,*]=coeffs[1]+(2*coeffs[2]*minutes)+(3*coeffs[3]*minutes^2)
  vel_funcs[i,*]=vel_funcs[i,*]*6.96e5/60  ;put in km/sec units
  acc_funcs[i,*]=2*coeffs[2]+(6*coeffs[3]*minutes)
  acc_funcs[i,*]=acc_funcs[i,*]*6.96e5/3600  ;put in km/sec^2 units
  h_tau=coeffs[0]+(coeffs[1]*t0)+(coeffs[2]*t0^2)+(coeffs[3]*t0^3)
  h_init=coeffs[0]+(coeffs[1]*xdata[0])+(coeffs[2]*xdata[0]^2)+ $
	(coeffs[3]*xdata[0]^3)
  h_final=coeffs[0]+(coeffs[1]*xdata[endpoint])+(coeffs[2]*xdata[endpoint]^2)+ $
	(coeffs[3]*xdata[endpoint]^3)
  minutes=minutes+offset

  ;
  ; Calculate initial, final, t0 velocities
  ;
  vel=coeffs(1)+(2*coeffs[2]*starttimes[i,*])+(3*coeffs[3]*starttimes[i,*]^2)

  ;initial velocity
  velinit=coeffs(1)+(2*coeffs[2]*xdata[0])+(3*coeffs[3]*xdata[0]^2) 
  vel=vel*6.96e5/60
  velinit=velinit*6.96e5/60
  velocities[i,0:1]=vel
  velocities[i,2]=velinit

  ; Final velocity
  vel=coeffs(1)+(2*coeffs[2]*xdata[endpoint])+(3*coeffs[3]*xdata[endpoint]^2)
  vel=vel*6.96e5/60
  velocities[i,3]=vel

  ; Velocity at t0
  vel=coeffs(1)+(2*coeffs[2]*t0)+(3*coeffs[3]*t0^2)
  vel=vel*6.96e5/60
  velocities[i,4]=vel

  ;
  ; Now it's time for the accelerations
  ;
  accel=((2*coeffs[2])+(6*coeffs[3]*(starttimes[i,*])))*6.96e5/3600
  accels[i,0:1]=accel
  initaccel=((2*coeffs[2])+(6*coeffs[3]*xdata[0]))*6.96e5/3600
  accels[i,2]=initaccel
  finalaccel=((2*coeffs[2])+(6*coeffs[3]*xdata[endpoint]))*6.96e5/3600
  accels[i,3]=finalaccel
  t0accel=((2*coeffs[2])+(6*coeffs[3]*t0))*6.96e5/3600

  PRINTF,unit,'==================== '+features[i]+' Data'+' ==================== ' 
  PRINTF,unit,''
  PRINTF,unit,'Returned coefficients: '
  PRINTF,unit,FORMAT='("Constant Term:",T16,e10.2," km or ",e10.2)',coeffs[0]*6.96e5,coeffs[0]
  PRINTF,unit,FORMAT='("Linear term:",T16,e10.2," km or ",e10.2)',coeffs[1]*6.96e5, coeffs[1]
  PRINTF,unit,FORMAT='("Quadratic term:",T16,e10.2," km or ",e10.2)',coeffs[2]*6.96e5, coeffs[2]
  PRINTF,unit,FORMAT='("Cubic term:",T16,e10.2," km or ",e10.2)',coeffs[3]*6.96e5,coeffs[3]
  PRINTF,unit,FORMAT='("Error in coeffs:",T16,e10.2," km or ",e10.2)',sigma*6.96e5,sigma
  PRINTF,UNIT,''
  PRINTF,UNIT,'********* Start time 1 **********'
  PRINTF,unit,FORMAT='("Start Time 1:",T26,A)',timestrarr[i,0]
  PRINTF,unit,FORMAT='("Start height:",T27,1F7.2)',startheights[i,0]
  PRINTF,unit,FORMAT='("Velocity:",T27,1F7.2)',velocities[i,0]
  PRINTF,unit,FORMAT='("Accel:",T22,1F12.4)',accels[i,0]
  PRINTF,unit,''
  PRINTF,unit,'********* Start time 2 **********'
  PRINTF,unit,FORMAT='("Start Time 2:",T26,A)',timestrarr[i,1]
  PRINTF,unit,FORMAT='("Start height:",T27,1F7.2)',startheights[i,1]
  PRINTF,unit,FORMAT='("Velocity:",T27,1F7.2)',velocities[i,1]
  PRINTF,unit,FORMAT='("Accel:",T22,1F12.4)',accels[i,1]
  PRINTF,unit,''
  PRINTF,unit,'****** At First Data Point ******'
  PRINTF,unit,FORMAT='("Time:",T26,A)',format_time(xdata[0]+offset)
  PRINTF,unit,FORMAT='("Height:",T27,1F7.2)',h_init
  PRINTF,unit,FORMAT='("Velocity:",T27,1F7.2)',velocities[i,2]
  PRINTF,unit,FORMAT='("Accel:",T22,1F12.4)',accels[i,2]
  PRINTF,unit,''
  PRINTF,unit,'******* At last datapoint *******'
  PRINTF,unit,FORMAT='("Time:",T26,A)',format_time(xdata[endpoint]+offset)
  PRINTF,unit,FORMAT='("Height:",T27,1F7.2)',h_final
  PRINTF,unit,FORMAT='("Velocity:",T27,1F7.2)',velocities[i,3]
  PRINTF,unit,FORMAT='("Accel:",T22,1F12.4)',accels[i,3]
  PRINTF,unit,''
  PRINTF,unit,'************ At tau *************'
  PRINTF,unit,FORMAT='("Time:",T26,A)',timestrarr[i,2]
  PRINTF,unit,FORMAT='("Height:",T27,1F7.2)',h_tau
  PRINTF,unit,FORMAT='("Velocity:",T27,1F7.2)',velocities[i,4]
  PRINTF,unit,FORMAT='("Accel:",T22,1F12.4)',t0accel
  PRINTF,unit,''
  PRINTF,unit,''

  coeff_arr[i,*]=coeffs
  ENDFOR  ;}
  CLOSE,unit
  FREE_LUN,unit

  sz=SIZE(xdata)

  ;
  ; Cubic Plot========================================================
  ;
  colors=['BLACK','BLUE','RED','GREEN','YELLOW','ORANGE','PURPLE','MAGENTA', $
	  'TURQUOISE','BROWN']
  FOR c=0,9 DO BEGIN
    linecolor,c,colors[c]
  ENDFOR

  ;orig=MIN(starttimes[where(starttimes NE -1)])
  IF (NOT KEYWORD_SET(origin)) THEN orig=1000 ELSE BEGIN
    num=FIX(STR_SEP(origin,":"))
    orig=(num(0)*60.)+num(1)
  ENDELSE

  IF (NOT KEYWORD_SET(end_x)) THEN end_x=MAX(mins+10) ELSE BEGIN
    x=FIX(STR_SEP(end_x,":"))
    end_x=(x(0)*60)+x(1)
    IF (end_x LT MAX(mins)) THEN end_x=end_x+1440
  ENDELSE

  SET_PLOT,'PS'
  device,/color,/portrait,/inches,xsize=8.0,ysize=11.0,$
        ;xoffset=0.5,yoffset=3.5,bits_per_pixel=8,$
        xoffset=0,yoffset=0,bits_per_pixel=8,$
        filename='/tmp/idlC.ps'
 
  ;PLOT,mins,radii,BACKGROUND=255,COLOR=255, $

  ;
  ;  Trajectory Plot
  ;
  IF (KEYWORD_SET(noannotate)) THEN BEGIN
  PLOT,findgen(10),findgen(10),BACKGROUND=255,COLOR=0, $
    XRANGE=[orig,end_x], $
    YRANGE=[1,MAX(radii)], $
    POSITION=[.18,.43,.88,.93], $
    XTITLE='Time (UT)',YTITLE='Height (Solar Radii)',PSYM=1, $
    CHARSIZE=1.5,XTICKFORMAT='xticks'
  ENDIF ELSE BEGIN
  PLOT,findgen(10),findgen(10),BACKGROUND=255,COLOR=0, $
    TITLE='CME Trajectory Plot!C'+yrddoy+' - Cubic Fit', $
    XRANGE=[orig,end_x], $
    YRANGE=[1,MAX(radii)], $
    POSITION=[.18,.43,.88,.93], $
    XTITLE='Time (UT)',YTITLE='Height (Solar Radii)',PSYM=1, $
    CHARSIZE=1.0,XTICKFORMAT='xticks'
  ENDELSE

  sym=[1,2,4,5,6]
  FOR p=0,num_features-1 DO BEGIN
    errors=SQRT(data[where(data.feature EQ features[p])].error)
    xdata=data[where(data.feature EQ features[p])].mins
    ydata=data[where(data.feature EQ features[p])].radius
    OPLOTERR,xdata,ydata,errors
    OPLOT,xdata,ydata,COLOR=p,PSYM=sym[p] ;,SYMSIZE=0.5
    ;xrng=WHERE(functions[p,*] GE MIN(startheights)) 
    xrng=WHERE(functions[p,*] GE 0.5) 
    OPLOT,minutes[xrng],functions[p,xrng],COLOR=p
    OPLOT,starttimes[p,*]+offset,startheights[p,*],COLOR=3,PSYM=5
  ENDFOR 

  IF (NOT(KEYWORD_SET(noannotate))) THEN BEGIN  ;{
  FOR d=0,num_datasets-1 DO BEGIN
    x_shift=0.2+(d/20.)
    XYOUTS,x_shift,0.9,datasets[d],/NORMAL
  ENDFOR

  FOR j=0, num_features-1 DO BEGIN  ;{
  y_shift=0.35-(j/11.)
  XYOUTS,0.5,y_shift,features[j]+' ('+colors[j]+')', $
	   /NORMAL,ALIGNMENT=0.5,COLOR=j

  XYOUTS,0.15,y_shift-.02,'Start time 1 (UT):', /NORMAL,ALIGNMENT=0.0, COLOR=j
  XYOUTS,0.45,y_shift-.02,timestrarr[j,0], $
	 /NORMAL,ALIGNMENT=1.0, COLOR=j
  XYOUTS,0.55,y_shift-.02,'Start time 2 (UT):', /NORMAL,ALIGNMENT=0.0, COLOR=j
  XYOUTS,0.9,y_shift-.02,timestrarr[j,1], $
	 /NORMAL,ALIGNMENT=1.0, COLOR=j

  XYOUTS,0.15,y_shift-.035,'Start Height (Rsun):',/NORMAL, $
	 ALIGNMENT=0.0,COLOR=j
  XYOUTS,0.45,y_shift-.035,STRING(startheights[j,0],FORMAT="(f6.1)"), $
	/NORMAL,ALIGNMENT=1.0,COLOR=j

  XYOUTS,0.55,y_shift-.035,'Start Height (Rsun):',/NORMAL, $
	 ALIGNMENT=0.0,COLOR=j
  XYOUTS,0.9,y_shift-.035,STRING(startheights[j,1],FORMAT="(f6.1)"), $
	/NORMAL,ALIGNMENT=1.0,COLOR=j


  XYOUTS,0.15,y_shift-.050,'Velocity (km/s):',/NORMAL, $
	 ALIGNMENT=0.0,COLOR=j
  XYOUTS,0.45,y_shift-.050,STRING(velocities[j,0],FORMAT="(f10.2)"), $
	/NORMAL,ALIGNMENT=1.0,COLOR=j
  XYOUTS,0.55,y_shift-.05,'Velocity (km/s):',/NORMAL, $
	ALIGNMENT=0.0,COLOR=j
  XYOUTS,0.9,y_shift-.05,STRING(velocities[j,1],FORMAT="(f10.2)"), $
	/NORMAL,ALIGNMENT=1.0,color=j

  XYOUTS,0.15,y_shift-.065,'Accel (km/s!U2!N):',/NORMAL, $
	ALIGNMENT=0.0,COLOR=j
  XYOUTS,0.45,y_shift-.065,STRING(accels[j,0],FORMAT="(f8.4)"), $
	/NORMAL,ALIGNMENT=1.0,COLOR=j
  XYOUTS,0.55,y_shift-.065,'Accel (km/s!U2!N):',/NORMAL, $
	ALIGNMENT=0.0,COLOR=j
  XYOUTS,0.9,y_shift-.065,STRING(accels[j,1],FORMAT="(f8.4)"), $
	/NORMAL,ALIGNMENT=1.0,COLOR=j

  ENDFOR ;}
  ENDIF  ;}

  ; 
  ;  Velocity Plot 
  ; 
  IF (KEYWORD_SET(noannotate)) THEN BEGIN
  PLOT,findgen(10),findgen(10),BACKGROUND=255,COLOR=0, $
    XRANGE=[orig,end_x], $
    YRANGE=[0,MAX(velocities)], $
    POSITION=[.18,.43,.88,.93], $
    XTITLE='Time (UT)',YTITLE='Velocity (km/sec)',PSYM=1, $
    CHARSIZE=1.5,XTICKFORMAT='xticks'
  ENDIF ELSE BEGIN
  PLOT,findgen(10),findgen(10),BACKGROUND=255,COLOR=0, $
    TITLE='CME Velocity Plot!C'+yrddoy+' - Cubic Fit', $
    XRANGE=[orig,end_x], $
    YRANGE=[0,MAX(velocities)], $
    POSITION=[.18,.43,.88,.93], $
    XTITLE='Time (UT)',YTITLE='Velocity (km/sec)',PSYM=1, $
    CHARSIZE=1.0,XTICKFORMAT='xticks'
  ENDELSE

  FOR p=0,num_features-1 DO BEGIN
    xrng=WHERE(vel_funcs[p,*] GE -100.5) 
    OPLOT,minutes[xrng],vel_funcs[p,xrng],COLOR=p 
    OPLOT,starttimes[p,*]+offset,startheights[p,*],COLOR=3,PSYM=5
  ENDFOR

  ;
  ;  Accel plots
  ;
  zero=minutes
  zero[*]=0
  xaxis=minutes+orig
  IF (KEYWORD_SET(noannotate)) THEN BEGIN
  PLOT,findgen(10),findgen(10),BACKGROUND=255,COLOR=0, $
    XRANGE=[orig,end_x], $
    YRANGE=[min(accels),MAX(accels)], $
    POSITION=[.18,.43,.88,.93], $
    XTITLE='Time (UT)',YTITLE='Acceleration (km/sec^2)',PSYM=1, $
    CHARSIZE=1.5,XTICKFORMAT='xticks'
  OPLOT,xaxis,zero,COLOR=0
  ENDIF ELSE BEGIN
  PLOT,findgen(10),findgen(10),BACKGROUND=255,COLOR=0, $
    TITLE='CME Acceleration Plot!C'+yrddoy+' - Cubic Fit', $
    XRANGE=[orig,end_x], $
    YRANGE=[MIN(accels),MAX(accels)], $
    POSITION=[.18,.43,.88,.93], $
    XTITLE='Time (UT)',YTITLE='Acceleration (km/sec)',PSYM=1, $
    CHARSIZE=1.0,XTICKFORMAT='xticks'
  OPLOT,xaxis,zero,COLOR=0
  ENDELSE
 
  FOR p=0,num_features-1 DO BEGIN
    xrng=WHERE(acc_funcs[p,*] GE -10.5)
    OPLOT,minutes[xrng],acc_funcs[p,xrng],COLOR=p
    OPLOT,starttimes[p,*]+offset,startheights[p,*],COLOR=3,PSYM=5
  ENDFOR 


  DEVICE,/CLOSE
  SET_PLOT,'X'

RETURN, coeffs
END