; ; This function performs a 1st 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_fit1,data 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+'.fit1.out' OPENW,unit,filename,/GET_LUN PRINTF,unit,'' PRINTF,unit,'Order Fit: 1' PRINTF,unit,'' ; ; Linear Fit ========================================================= ; coeff_arr=FLTARR(num_features,2) starttimes=FLTARR(num_features) velocities=FLTARR(num_features) ;functions=FLTARR(num_features,num_meas) 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,1,yfit,yband,sigma,corr) starttime=(1.0-coeffs(0))/coeffs(1) ;1st order start time (mins) time_str=format_time(starttime) velocity=coeffs(1)*(6.96e5)/60 ;Radial velocity ;functions[i,*]=coeffs(0)+(coeffs(1)*mins) functions[i,*]=coeffs(0)+(coeffs(1)*minutes) starttimes[i]=starttime velocities[i]=velocity 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,'Error in coeffs: '+STRING(sigma*6.96e5)+' km OR ',sigma PRINTF,UNIT,'' PRINTF,unit,'Start Time: '+time_str PRINTF,unit,'Radial Velocity: '+STRING(velocity)+' km/sec' PRINTF,unit,'' print,features[i]+' start time: '+time_str coeff_arr[i,*]=coeffs ENDFOR ;} CLOSE,unit FREE_LUN,unit sz=SIZE(xdata) ; ; Linear Plot======================================================== ; colors=['BLACK','BLUE','RED','GREEN','YELLOW','ORANGE','PURPLE','MAGENTA', $ 'TURQUOISE','BROWN'] FOR c=0,num_features-1 DO BEGIN linecolor,c,colors[c] ENDFOR orig=MIN(starttimes) end_x=MAX(mins+10) 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, $ PLOT,findgen(10),findgen(10),BACKGROUND=255,COLOR=0, $ TITLE='CME Trajectory Plot!C'+yrddoy+' - Linear Fit', $ XRANGE=[orig,end_x], $ YRANGE=[1,MAX(radii)], $ POSITION=[.18,.45,.88,.95], $ XTITLE='Time (UT)',YTITLE='Height (Solar Radii)',PSYM=1, $ CHARSIZE=1.0,XTICKFORMAT='xticks' 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 1.0) OPLOT,minutes[xrng],functions[p,xrng],COLOR=p 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.2,y_shift-.02,'Start time (UT):', /NORMAL,ALIGNMENT=0.0, COLOR=j XYOUTS,0.7,y_shift-.02,format_time(starttimes[j]), $ /NORMAL,ALIGNMENT=1.0, COLOR=j XYOUTS,0.2,y_shift-.04,'Average Velocity (km/sec):',/NORMAL, $ ALIGNMENT=0.0,COLOR=j XYOUTS,0.7,y_shift-.04,velocities[j],/NORMAL,ALIGNMENT=1.0,COLOR=j ENDFOR DEVICE,/CLOSE SET_PLOT,'X' RETURN, coeffs END