;
;  This function performs a 2nd 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_fit2,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/'
;  OPENW,unit,path+filename,/GET_LUN
;  instead of hardwiring the path use current working directory: J.Burkepile June 2008

  filename=yrddoy+'.fit2.out'
  OPENW,unit,filename,/GET_LUN
  PRINTF,unit,''
  PRINTF,unit,'Order Fit:  2'
  PRINTF,unit,''

  ;
  ; Quadratic Fit =========================================================
  ;
  coeff_arr=FLTARR(num_features,3)
  starttimes=FLTARR(num_features,1) 
  timestrarr=STRARR(num_features,1)
  startheights=FLTARR(num_features,1)
  velocities=FLTARR(num_features,2)  ;initial, final velocity for each feat.
  accels=FLTARR(num_features,1)
  functions=FLTARR(num_features,sz(1))

  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=data[where(data.feature EQ features[i])].error

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

  ;
  ; Calculate start times
  ;  If quadratic solition in imaginary, then parabola does not intersect
  ;  radius=1, set derivative = 1 and find minimum.

  IF (coeffs[1]^2 LT (4*coeffs[2]*(coeffs[0]-1.0))) THEN BEGIN
    starttimes[i,0]=-coeffs[1]/(2*coeffs[2])
  ENDIF ELSE BEGIN

  ; If quadratic is real, then parabola intersects radius=1 and this is
  ;  our start time for the event.
    quad1=-coeffs[1]/(2*coeffs[2])
    quad2=SQRT(coeffs[1]^2-4*coeffs[2]*(coeffs[0]-1.0))/(2*coeffs[2])
    starttimes[i,0]=quad1+quad2
  ENDELSE

  timestrarr[i,0]=format_time(starttimes[i,0])

  ;
  ;Start height is the quadratic evaluated at the start time.
  ;
  startheights[i,0]=coeffs(0)+(coeffs(1)*starttimes[i,0])+ $
	(coeffs(2)*starttimes[i,0]^2)
  index=where(starttimes EQ -1)
  IF (index[0] NE -1) THEN startheights[index] = 0.0
  functions[i,*]=coeffs(0)+(coeffs(1)*minutes)+coeffs(2)*minutes^2

  ;
  ; Calculate initial and final velocities
  ;
  vel1=coeffs(1)+(2*coeffs[2]*starttimes[i,0])
  vel1=vel1*6.96e5/60
  vel2=coeffs(1)+(2*coeffs[2]*MAX(mins))
  vel2=vel2*6.96e5/60
  velocities[i,0]=vel1
  velocities[i,1]=vel2

  ;
  ; Now it's time for the accelerations
  ;
  accel=2*coeffs[2]*6.96e5/3600
  accels[i,0]=accel

  PRINTF,unit,features[i]+' Data'
  PRINTF,unit,'=================================================='
  PRINTF,unit,''
  PRINTF,unit,'Returned coefficients: '
  PRINTF,unit,'Constant term:  '+STRING(coeffs[0]*6.96e5)+' km OR ',coeffs[0]
  PRINTF,unit,'Linear term:    '+STRING(coeffs[1]*6.96e5)+' km/sec OR ', $
	coeffs[1]
  PRINTF,unit,'Quadratic term: '+STRING(coeffs[2]*6.96e5)+' km/sec^2 OR ', $
	coeffs[2]
  PRINTF,unit,'Error in coeffs: '+STRING(sigma*6.96e5)+' km OR ',sigma
  PRINTF,UNIT,''
  PRINTF,unit,'Start Time:  '+timestrarr[i,0]
  PRINTF,unit,''
  PRINTF,unit,'Start height: '+STRING(startheights[i,0])
  PRINTF,unit,''
  PRINTF,unit,'Initial Velocity: '+STRING(velocities[i,0])
  PRINTF,unit,'Final Velocity: '+STRING(velocities[i,1])
  PRINTF,unit,''
  PRINTF,unit,'Average Accel: '+STRING(accels[i,0])

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

  sz=SIZE(xdata)

  ;
  ; Quadratic 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=1100 ELSE BEGIN
  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, $
  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+' - Quadratic 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
    xdata=data[where(data.feature EQ features[p])].mins
    ydata=data[where(data.feature EQ features[p])].radius
    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,*],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
  ENDIF

  IF (NOT(KEYWORD_SET(noannotate))) THEN BEGIN  ;{
  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 (UT):', /NORMAL,ALIGNMENT=0.0, COLOR=j
  XYOUTS,0.45,y_shift-.02,format_time(starttimes[j,0]), $
	 /NORMAL,ALIGNMENT=1.0, COLOR=j

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

  XYOUTS,0.15,y_shift-.050,'Initial 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.15,y_shift-.065,'Final Velocity (km/s):',/NORMAL, $
	ALIGNMENT=0.0,COLOR=j
  XYOUTS,0.45,y_shift-.065,STRING(velocities[j,1],FORMAT="(f10.2)"), $
	/NORMAL,ALIGNMENT=1.0,color=j

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

  ENDFOR
  ENDIF  ;}

  DEVICE,/CLOSE
  SET_PLOT,'X'

RETURN, coeffs
END