ridge/scripts/polyfit.pro
Roman Krivonos 2552bdb552 minor
2024-07-06 17:35:25 +03:00

59 lines
1.3 KiB
Prolog

pro polyfit
green = GETCOLOR('green', 100)
red = GETCOLOR('red', 100)
enkey="B21"
fn = '../products/detcnts.'+enkey+'.ALL.resid.fits'
table = readfits(fn, header, EXTEN_NO=1, /SILENT)
;;tstart = TBGET( header, table, 'TSTART', /NOSCALE)
;;tstop = TBGET( header, table, 'TSTOP', /NOSCALE)
rev = TBGET( header, table, 'REV', /NOSCALE)
data = TBGET( header, table, 'RESID', /NOSCALE)
time = LINDGEN(N_ELEMENTS(data))
rev_min=min(rev)
rev_max=max(rev)
plot_x=[]
plot_y=[]
for r=rev_min, rev_max do begin
index = where(rev eq r, count)
if(count lt 10) then continue
mn = AVG(data(index))
push,plot_x,r
push,plot_y,mn
print,r,count,mn
endfor
plot,plot_x,plot_y ;, yrange=[-0.0001,0.0001]
return
;;COEFF = ROBUST_LINEFIT(time, rate, YFIT, SIG, COEF_SIG)
COEFF = ROBUST_POLY_FIT(time,data,7,YFIT,SIG, /DOUBLE)
print,SIG
print,COEFF
;mean_polyfit=AVG(YFIT)
;print,'mean fit',mean_polyfit
;print,'resid',SIG/mean_polyfit*100
poly = coeff;/mean_polyfit
;;save,FILENAME = enkey+'_poly.sav', poly, mean_polyfit
scale=1.2
;;yfit0 = (poly[0] + poly[1]*time + poly[2]*time^2 + poly[3]*time^3)*scale
plot,time,data;, yrange=[-0.0001,0.0001]
oplot,time,YFIT,color=green
;;oplot,time,YFIT0,color=red
end