-
Notifications
You must be signed in to change notification settings - Fork 0
/
create_traectory.pro
98 lines (86 loc) · 2.85 KB
/
create_traectory.pro
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
function create_traectory,ima,NP=Np,X_BEG=x_beg,X_END=X_END,WX=wx,NDEG=Ndeg,PLOT=plot
;ïîñòðîåíèå òðàåêòîðèè äëÿ èçîáðàæåíèÿ FLAT ñ ìàñêîé SCORPIO-2
;Np - ÷èñëî òðàåêòîðèé
;X_beg - êîîðäèíàòà ïåðâîãî ñòðîáà âäîëü ùåëè
;WX - øèðèíà ñòðîáà
;Ndeg - ñòåïåíü ïîëèíîìà àïïðîêñèìàöèè òðàåêòîðèè
!P.multi=[0,1,1]
a=size(ima) & Nx=a(1) & Ny=a(2)
if not(keyword_set(Np)) then Np=10
if not(keyword_set(X_beg)) then X_beg=wx
if not(keyword_set(X_end)) then X_end=Nx-1
if not(keyword_set(WX)) then wx=80
if not(keyword_set(Ndeg)) then Ndeg=3
Npos=fix((x_end-x_beg )/wx/2)
xpos=findgen(Npos)*2*wx+x_beg
print,xpos
;print,Npos,Ny
;print,xpos
y=findgen(Ny)
vector_y=fltarr(Npos,Ny)
ypk=fltarr(Np,Npos)
;ôîðìèðîâàíèå ðàçðåçîâ (ñòðîáîâ)
for k=0,Npos-1 do vector_y(k,*)=total(ima(xpos(k)-wx:xpos(k)+wx,*),1)
;pk=find_peaks(vector_y(0,*),W=10,/plot)
;return,0
;îïðåäåëåíèå êîîðäèíàò ïèêîâ ñðåäíåãî ðàçðåçà
mpk=find_peaks(vector_y(Npos/2,*),W=5,tresh=1000,/plot)
print,mpk
if N_elements(mpk) ne Np then begin
print,'number peaks not equal '+string(Np,format='(I2)')+'!!!'
RETURN,0
endif
;ypk(*,Npos/2)=pk
print,' mean',mpk
;îïðåäåëåíèå ïîçèöèé âñåõ ïèêîâ
for k=0,Npos-1 do begin
pk=find_peaks(vector_y(k,*),W=5,tresh=1000,/plot)
print,k,N_elements(pk)
;print,pk
index=compare_list(mpk,pk,20)
ypk(index(*,1),k)=pk(index(*,1))
endfor
;ïîñòðîåíèå òðàåêòîðèé flat
x=findgen(Nx)
tra=fltarr(Nx,Np)
wdelete,2
if keyword_set(plot) then begin
window,22,xpos=0,ypos=0,title='SPECTRA TRAECTORY'
;set_plot,'PS'
;device,file='g:\traectory.eps',/portrait,/encapsulated,xsize=16,ysize=10
plot,[0,Nx],[0,Ny],/nodata,xst=1,yst=1,title='Image binning 2x2 px',$
;plot,[0,Nx],[50,450],/nodata,xst=1,yst=1,yminor=5,title='Image binning 2x2 px',$
ytitle='Distance along slit, px',xtitle='Distance along dispersion, px',charsize=1
endif
for j=0,Npos-1 do begin
R=where(ypk(*,j) ne 0, ind)
if ind ne 0 and ind ne Np then ypk(*,j)=0
endfor
for k=0,Np-1 do begin
;remove zero values
R=where(ypk(k,*) ne 0,ind)
print,R
f=goodpoly(xpos(R),ypk(k,R),Ndeg,3,Yfit)
robomean,ypk(k,R)-Yfit,3,0.5,avg_fit,rms_fit
for j=0,Ndeg do tra(*,k)=tra(*,k)+f(j)*x^j
if keyword_set(plot) then begin
oplot,xpos(R),ypk(k,R),psym=6,symsize=0.6
oplot,x,tra(*,k)
xyouts,Nx/2,tra(Nx/2,k)+10,'rms ='+string(rms_fit,format='(F5.2)')+' px',/data,align=0.5
endif
endfor
;device,/close
;set_plot,'WIN'
return,tra
end
;LOGFILE=DIALOG_PICKFILE(/READ, FILTER = '*.txt',path='h:\spectraPOL.log')
;wdir=sxpar(read_table(LOGFILE),'w_dir')
path='h:\red_data.pol\AGN\'
LOGFILE=DIALOG_PICKFILE(/READ, FILTER = '*.txt',path=path+'LOGS')
wdir=str_sep(sxpar(read_table(LOGFILE),'w_dir'),'\')
wdir=path+wdir(N_elements(wdir)-2)+'\'
flat=readfits(wdir+'eta_i.fts')
R=where(flat lt 0) & flat(R)=0
tra=create_traectory(flat,NP=3,WX=20,NDEG=2,X_beg=200,/plot)
;writefits,wdir+'traectory.fit',tra
end