Документ взят из кэша поисковой машины. Адрес оригинального документа : http://lnfm1.sai.msu.ru/~kupry/SOFTWARE/mfsff_v12.3_winyuk.pro Дата изменения: Sun Oct 9 23:32:42 2011 Дата индексирования: Mon Oct 1 23:03:34 2012 Кодировка: Поисковые слова: п п п п р п р п р п р п

;
;
; ++++++++++++++++++++++++++++++++++++++++++++++++++++++++
; + +
; + FUNCTIONS & PROCEDURES +
; + +
; ++++++++++++++++++++++++++++++++++++++++++++++++++++++++
;
;
pro make_setup
;
common setup_cb,idl_4degpolyfitting_sw,polyfitting_path
;
; procedure doing only a set-up used throughout whole code
; set-up variable are tresnfered into procedures and main
; code via common block setup_cb
;
; if idl_4degpolyfitting_sw equal to zero, code written
; and compiled in FORTRAN is used for fitting with 3th degree
; polynome for estimation of the spectra curvature. If it is
; greater than zero, idl procedures from the file polyfit4deg.pro
; are used.
; The IDL procedure has lower numerical stability than the FORTRAN code!!!
idl_4degpolyfitting_sw=1
;
; path the compiled FORTRAN code (idl_4degpolyfitting_sw=0) or to
; directory containing the IDL procedure polyfit4deg.pro (idl_4degpolyfitting_sw=1).
; It is highly recommended that the FORTRAN code is compiled using
; compile script under UNIX systems
;
polyfitting_path='c:\MFS_zpracovani\software' ; <--- for IDL procedure
; polyfitting_path='./polyfit.out' ; <--- for FORTRAN CODE
;
;
;
;if(idl_4degpolyfitting_sw gt 0) then begin
; !PATH=!PATH+':'+polyfitting_path
;endif
;
end
;
;
;
function poly4fnc,x,a
;
y=a[0]+a[1]*x+a[2]*x^2+a[3]*x^3
yder0=dblarr(n_elements(x))+1.0d0
yder1=x
yder2=x^2
yder3=x^3
return,[ [y], [yder0], [yder1], [yder2] ,[yder3] ]
end
;
;
;
pro polyfit4deg,x,y,params,meas_err=meas_err,winn=winn
;
if(!version.os_family eq 'unix') then windev='X' else windev='win'
origdev=!d.name
tvlct,rcol_orig,gcol_orig,bcol_orig,/get
set_plot,windev
;
rcol=[0.,1.,0.,1.]*255
gcol=[0.,1.,1.,0.]*255
bcol=[0.,1.,0.,0.]*255
tvlct,rcol,gcol,bcol
;
device,dec=0
if(n_elements(meas_err) eq 0) then meas_err=sqrt(y)
if(n_elements(winn) gt 0) then begin
window,winn
endif else begin
window,/free
winn=!window
endelse
a=[1.d0,1.d0,1.d0,1.d0]
ftng_agn_ln:
plot,x,y,/nodata
oplot,x,y,psym=1,color=3
tmpstr=' '
print,''
print,''
print,'function: a0+a1*x+a2*x^2+a3*x^3'
print,'a0=',a[0],'(leave blank for no change)'
read,tmpstr
if(strlen(tmpstr) gt 0) then a[0]=double(tmpstr)
print,''
print,'a1=',a[1],'(leave blank for no change)'
read,tmpstr
if(strlen(tmpstr) gt 0) then a[1]=double(tmpstr)
print,' '
print,'a2=',a[2],'(leave blank for no change)'
read,tmpstr
if(strlen(tmpstr) gt 0) then a[2]=double(tmpstr)
print,' '
print,'a3=',a[3],'(leave blank for no change)'
read,tmpstr
if(strlen(tmpstr) gt 0) then a[3]=double(tmpstr)
print,' '
print,' '
print,'Choose whether a parameter is taken in the fitting'
print,' Otherwise it is kept constant'
print,' '
pmtft=[0,0,0,0]
tmpstr=' '
print,'Take a[0] into fitting? (Y/n)'
read,tmpstr
tmpstr=strupcase(tmpstr)
if(tmpstr ne 'N') then pmtft[0]=1
print,' '
tmpstr=' '
print,'Take a[1] into fitting? (Y/n)'
read,tmpstr
tmpstr=strupcase(tmpstr)
if(tmpstr ne 'N') then pmtft[1]=1
print,' '
tmpstr=' '
print,'Take a[2] into fitting? (Y/n)'
read,tmpstr
tmpstr=strupcase(tmpstr)
if(tmpstr ne 'N') then pmtft[2]=1
print,' '
tmpstr=' '
print,'Take a[3] into fitting? (Y/n)'
read,tmpstr
tmpstr=strupcase(tmpstr)
if(tmpstr ne 'N') then pmtft[3]=1
print,' '
fitxvec=lmfit(x,y,a,chisq=ch_sq,/double,fita=pmtft,measure_errors=meas_err,$
function_name='poly4fnc',sigma=sgm_a)
print,'function: a0+a1*x+a2*x^2+a3*x^3'
for i=0,3 do begin
aistr='a['+strcompress(string(i),/remove_all)+']='+strcompress(string(a[i]),/remove_all)
if(pmtft[i] eq 1) then print,aistr,' +/-',abs(sgm_a[i]/a[i]*100.),' %' $
else print,aistr,' +/- 0.0 %'
endfor
xxstep=(max(x)-min(x))/100.
xxx=findgen(101)*xxstep+min(x)
fitvec=a[0]+a[1]*xxx+a[2]*xxx^2+a[3]*xxx^3
window,winn
plot,x,y,/nodata
oplot,x,y,psym=1,color=3
oplot,xxx,fitvec,color=2
print,''
print,'--------------------------------'
print,'Continue with fitting? (Y/n)'
print,' '
sw_cft=' '
read,sw_cft
sw_cft=strupcase(sw_cft)
if(sw_cft ne 'N') then goto,ftng_agn_ln
;
params=a
tvlct,rcol_orig,gcol_orig,bcol_orig
set_plot,origdev
end
;
;
;
function dtfilter,mat,threshold=threshold
if(n_elements(threshold) eq 0) then threshold=2
if(threshold lt 1.) then begin
resmat=mat
endif else begin
top_limit=double(threshold)
bottom_limit=1.d0/double(threshold)
resmat=((double(mat)bottom_limit)
endelse
return,resmat
end
;
;
;
pro rmvar,var
var1=temporary(var)
end
;
;
;
function ysection,y1,y2,step
y1s=min([y1,y2])
y2s=max([y1,y2])
rn_resvec=(y2s-y1s)/step+1.
n_resvec=fix(rn_resvec)
if((rn_resvec-float(n_resvec)) gt 0.) then n_resvec1=n_resvec+1L else $
n_resvec1=n_resvec
resvec=fltarr(n_resvec1)
resvec[0:(n_resvec-1L)]=findgen(n_resvec)*step+y1s
if(n_resvec1 gt n_resvec) then resvec[n_resvec]=y2
;
return,resvec
end
;
;
;
pro showimg,data=data
tvscl,data
end
;
;
;
pro sep_pathflnm,path_name,path,flnm,flext
; separate path, filename and extension from the whole path to file path_name
if(!version.os_family eq 'unix') then dir_separator='/'
if(!version.os_family eq 'Windows') then dir_separator='\'
posls=strpos(path_name,dir_separator,/reverse_search)
path=strmid(path_name,0,posls+1)
flnmwh=strmid(path_name,posls+1)
pospt=strpos(flnmwh,'.',/reverse_search)
flnm=strmid(flnmwh,0,pospt)
flext=strmid(flnmwh,pospt+1)
end
;
;
;
function hist_filter,inpmat,percentage=percentage
;
if(n_elements(percentage) eq 0) then percentage=10.d0
;
perc=float(percentage)/100.d0
h=histogram(inpmat,nbins=1000,locations=xh)
h=double(h)
xh=double(xh)
maxh=max(h)
lowlimhist=maxh*perc
indcs_abovelim=where(h ge lowlimhist)
if(indcs_abovelim[0] ne (-1)) then begin
ind_ll=min(indcs_abovelim)
ind_ul=max(indcs_abovelim)
ll=xh[ind_ll]
ul=xh[ind_ul]
outmat=((inpmat>ll) endif else begin
outmat=inpmat
endelse
;
return,outmat
end
;
;
;
pro lin_interp,xu,yu,x1,y1,extrapolate=extrapolate,outside_const=outside_const
;
;
;
if(not (keyword_set(extrapolate))) then extrapolate=0
if(not (keyword_set(outside_const))) then outside_const=0
if((extrapolate gt 0) and (outside_const gt 0)) then begin
print,'% lin_interp: Keywords extrapolate and outside_const'
print,'% cannot be used at the same time. Igno-'
print,'% ring both of them'
extrapolate=0 & outside_const=0
endif
;
;
nnn=n_elements(xu)
nny=n_elements(yu)
nnn1=n_elements(x1)
if (nnn ne nny) then begin
message,'the first two arguments must have the same dimensions!!! Exiting ...'
endif
if(nnn1 eq 0) then begin
message,'the third argument must have non-zero dimension!!! Exiting ...'
endif
;
;
x=xu
y=yu
y1=dblarr(nnn1)
;
if (nnn eq 1) then begin
for i=0,nnn1-1L do begin
if (x1[i] eq x[0]) then y1[i]=y[0]
endfor
goto,lenl
endif
;
; sorting
sortagl:nswaps=0L
for i=0,nnn-2L do begin
if(x[i] gt x[i+1]) then begin
xtmp=x[i]
x[i]=x[i+1]
x[i+1]=xtmp
ytmp=y[i]
y[i]=y[i+1]
y[i+1]=ytmp
nswaps=nswaps+1L
endif
endfor
if (nswaps gt 0) then goto,sortagl
;
;
; linear interpolation point-by-point
;
j=0
for i=0,nnn1-1L do begin
k_coef=0.d0
b_coef=0.d0
case 1 of
(x1[i] lt x[0]):begin
if(extrapolate gt 0) then begin
k_coef=(y[1]-y[0])/(x[1]-x[0])
b_coef=y[1]-k_coef*x[1]
endif
if(outside_const gt 0) then begin
k_coef=0.d0
b_coef=y[0]
endif
end
(x1[i] gt x[nnn-1L]):begin
if(extrapolate gt 1) then begin
k_coef=(y[nnn-1L]-y[nnn-2L])/(x[nnn-1L]-x[nnn-2L])
b_coef=y[nnn-1L]-k_coef*x[nnn-1L]
endif
if(outside_const gt 0) then begin
k_coef=0.d0
b_coef=y[nnn-1L]
endif
end
else:begin
znzn:
if (x1[i] eq x[j]) then begin
y1[i]=y[j]
goto,eorcl
endif
if (x1[i] lt x[j]) then begin
j=j-1L
goto,znzn
endif
if (x1[i] gt x[j+1L]) then begin
j=j+1L
goto,znzn
endif
if (x1[i] eq x[j+1L]) then begin
y1[i]=y[j+1L]
goto,eorcl
endif
k_coef=(y[j+1]-y[j])/(x[j+1]-x[j])
b_coef=y[j+1L]-k_coef*x[j+1L]
end
endcase
y1[i]=k_coef*x1[i]+b_coef
eorcl:
endfor
;
;
lenl:
end
;
;
;
pro corr_incl,mat1,mat2,tanhorizincl
;
; correction for vertical inclination of the spectrum
;
; mat1 .................... input matrix
; mat2 .................... output matrix
; tanhorizincl ............ tangens of inclination angle measured
; from X-axis in counter-clockwise direction
;
xdim=n_elements(mat1[*,0])
ydim=n_elements(mat1[0,*])
mat2=dblarr(xdim,ydim)*0.0d0
ypositions=findgen(ydim)
;
for ix=0,xdim-1L do begin
yshift=double(ix)*tanhorizincl
yvalvec=double(reform(mat1[ix,*]))
lin_interp,ypositions-yshift,yvalvec,ypositions,yvalvec1
mat2[ix,*]=yvalvec1
endfor
mm=[min(mat1),max(mat1)]
indcsorng=where(mat2 lt mm[0])
if(indcsorng[0] ne (-1)) then mat2[indcsorng]=mm[1]
mat2=(mat2 end
;
;
;
pro corr_curv,mat1,mat2,polyparams
;
; correction for curvature of the spectra estimated
; from horizontal curvature of chosen spectral line
;
; mat1 .................... input matrix
; mat2 .................... output matrix
; polyparams .............. vector of coefficient of
; polynomial fnc sum_i=0^N a(i)*x^i,
; expressing the curvature
;
xdim=n_elements(mat1[*,0])
ydim=n_elements(mat1[0,*])
mat2=dblarr(xdim,ydim)*0.0d0
xpositions=findgen(xdim)
n_polyparams=n_elements(polyparams)
for iy=0,ydim-1L do begin
xshift=0.d0
for j=0,n_polyparams-1L do begin
if(j eq 0) then val1=polyparams[0] else val1=polyparams[j]*double(iy)^double(j)
xshift=xshift+val1
endfor
xvalvec=double(reform(mat1[*,iy]))
lin_interp,xpositions-xshift,xvalvec,xpositions,xvalvec1
mat2[*,iy]=xvalvec1
endfor
mm=[min(mat1),max(mat1)]
indcsorng=where(mat2 lt mm[0])
if(indcsorng[0] ne (-1)) then mat2[indcsorng]=mm[1]
mat2=(mat2 end
;
;
;
function numb2string,numb
str1=string(numb)
return,strcompress(str1,/remove_all)
end
;
;
;
pro estim_ffshift,mat,ffmat,shift,shift_range=shift_range,shift_step=shift_step,$
correlation=correlation,auto=auto,pltwindow=pltwindow
;
; DESCRIPTION:
; estimate shift of the slit-flat frame according to
; spectrum frame. For estimation of the best agreement
; between shifted slit-flat and spectrum frames the
; cross-correlation is used
;
; INPUT:
; mat ............ observed spectr-um/-a
; ffmat .......... slit-flat frame
;
; OUTPUT:
; shift .......... estimated y-shift in pxl, that the slit-flat
; frame should shifted by.
;
; KEYWORDS
; shift_rangea ..... two-member vector containg range of shift value,
; in which the best shift is searched for. If this
; keyword is omitted, default range +/- 10 pxls is
; taken.
;
;
; Written by P. Schwartz, July the 4, 2011
;
;++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
;
;
;
if(n_elements(auto) eq 0) then auto=0
if(n_elements(shift_range) lt 2) then shrng=[-10.,10.] else $
shrng=[min(shift_range),max(shift_range)]
if(n_elements(shift_step) eq 0) then begin
shift_step=1.
endif else begin
if(shift_step eq 0.) then shift_step=1. else shift_step=abs(shift_step)
endelse
if(n_elements(pltwindow) eq 0) then pltwindow=0
;
dimensions=size(mat)
xdim=dimensions[1]
ydim=dimensions[2]
ypositions=findgen(ydim)
if (dimensions[0] eq 3) then begin
zdim=dimensions[3]
avgmat=dblarr(xdim,ydim)*0.d0
for iz=0L,zdim-1L do avgmat=avgmat+reform(mat[*,*,iz])
avgmat=avgmat/double(zdim)
endif else begin
avgmat=mat
endelse
;
imat_yav=0L
iffmat_yav=0L
mat_yav=dblarr(ydim)*0.d0
ffmat_yav=dblarr(ydim)*0.d0
for ix=0L,xdim-1L do begin
tmpvec=reform(mat[ix,*])
ii_zer=total(abs(tmpvec))
if(ii_zer ne 0.) then begin
mat_yav=mat_yav+tmpvec
imat_yav=imat_yav+1L
endif
tmpvec=reform(ffmat[ix,*])
ii_zer=total(abs(tmpvec))
if(ii_zer ne 0.) then begin
ffmat_yav=ffmat_yav+tmpvec
iffmat_yav=iffmat_yav+1L
endif
endfor
mat_yav=mat_yav/double(imat_yav)
tmpvec1=smooth(mat_yav,[11],/edge_truncate)
mat_yav1=mat_yav/tmpvec1
mat_yav_nm=mat_yav1-total(mat_yav1)/double(ydim)
mat_yav_nm=mat_yav_nm/(max(mat_yav_nm)-min(mat_yav_nm))
;
ffmat_yav=ffmat_yav/double(iffmat_yav)
ffmat_yav_nm=ffmat_yav-total(ffmat_yav)/double(ydim)
ffmat_yav_nm=ffmat_yav_nm/(max(ffmat_yav_nm)-min(ffmat_yav_nm))
;
if(auto gt 0) then begin
; vector of shifts
shftvec=findgen(long((shrng[1]-shrng[0])/shift_step)+1L)*shift_step+shrng[0]
n_shftvec=n_elements(shftvec)
if(shftvec[n_shftvec-1L] lt shrng[1]) then begin
n_shftvec=n_shftvec+1L
tmpvec2=dblarr(n_shftvec)
tmpvec2[0:n_shftvec-2L]=shftvec
tmpvec2[n_shftvec-1L]=shrng[1]
shftvec=temporary(tmpvec2)
endif
exzrinvec=where(shftvec eq 0.)
if(exzrinvec[0] eq (-1)) then begin
n_shftvec=n_shftvec+1L
tmpvec2=dblarr(n_shftvec)
indcs_ltzer=where(shftvec lt 0.)
indcs_gtzer=where(shftvec gt 0.)
i4zer=0
if(indcs_ltzer[0] ne (-1)) then begin
nlt=n_elements(indcs_ltzer)
tmpvec2[0:nlt-1L]=shftvec[indcs_ltzer]
i4zer=nlt
endif
tmpvec2[i4zer]=0.0d0
if(indcs_gtzer[0] ne (-1)) then begin
ngt=n_elements(indcs_gtzer)
tmpvec2[i4zer+1L:*]=shftvec[indcs_gtzer]
endif
shftvec=temporary(tmpvec2)
endif
;
correlvec=dblarr(n_shftvec)
for isft=0L,n_shftvec-1L do begin
sftval=shftvec[isft]
lin_interp,ypositions,ffmat_yav_nm,ypositions+sftval,tmpvec3
correlvec[isft]=correlate(mat_yav_nm,tmpvec3,/double)
endfor
;
correlation=max(correlvec,ishmx)
shift=shftvec[ishmx]
lin_interp,ypositions,ffmat_yav_nm,ypositions+shift,tmpvec3
window,pltwindow
plot,mat_yav_nm
oplot,tmpvec3,line=1
endif else begin
print,' '
window,pltwindow
sftval=shift
ssall0546:
lin_interp,ypositions,ffmat_yav_nm,ypositions+sftval,tmpvec3
plot,mat_yav_nm
oplot,tmpvec3,line=1
correlval=correlate(mat_yav_nm,tmpvec3,/double)
print,'shift: ',sftval,' corelation: ',correlval
answ_s=' '
print,'Are you satisfied? (y/N)'
read,answ_s
answ_s=strupcase(answ_s)
if(answ_s ne 'Y') then begin
print,' '
print,'Change shift: '
read,sftval
print,' '
goto,ssall0546
endif
print,''
correlation=correlval
shift=sftval
print,' '
endelse
;
;
end
;
;
;
pro davidhalpprof,mu,wvl,int,relwvl=relwvl
;
;
mu=((mu<1.)>0.141)
centwvl=0.
imax=0.
nnn=0
;values of mu for 7 columns of intensities
musforcols=[1.0,0.8,0.6,0.436,0.312,0.28,0.141]
;
openr,1,'c:\MFS_zpracovani\software\david_halp.tab'
readf,1,centwvl,imax,nnn
; darkening coefficients
darcoef=fltarr(7)
intens=fltarr(7,nnn)
readf,1,darcoef
data=fltarr(8,nnn)
readf,1,data
close,1
for icol=1,7 do begin
intens(icol-1,*)=data(icol,*)/100.*imax*darcoef[icol-1]
endfor
halfprof=fltarr(2,nnn)
halfprof(0,*)=data(0,*)
complprof=fltarr(2,nnn*2-1)
for i=0,nnn-1 do begin
lin_interp,musforcols,reform(intens(*,i)),[mu],profval
halfprof(1,i)=profval
endfor
complprof(0,0:nnn-1)=reverse(halfprof(0,*),2)*(-1.0)
complprof(1,0:nnn-1)=reverse(halfprof(1,*),2)
complprof(*,nnn:nnn*2-2)=halfprof(*,1:nnn-1)
complprof(0,*)=complprof(0,*)+centwvl
lin_interp,complprof(0,*),complprof(1,*),wvl,int,/outside_const
relwvl=wvl-centwvl
;
end
;
;
;
pro find_circcent,points,center,radius=radius
;
; DESCRIPTION:
; Finds center of the circle defined passing through several points.
; Number of points should be even and larger than or equal than four.
; The more points is used better the center estimation is.
;
; INPUT:
; points ...... matrix 2xN, where N is number of points. It defines
; [x,y] coordinates of points which the circle is
; passing through. N must be larger than or equal to 4.
;
; OUTPUT:
; center ...... [x,y] vector of coordinates of the circle center
;
; OTIONAL OUTPUT:
; radius ...... keyword, is present and equal to some named variable,
; the estimated radius of the circle is stored into the
; variable
;
; Version 1.0 -- July the 11, 2011 by Pavol Schwartz
; schwartz@asu.cas.cz; pschwartz@ta3.sk
;
;=========================================================================
;
points=double(points)
npoints=n_elements(points[0,*])
center=[0.,0.]
if(npoints lt 4) then begin
print,'% find_circcent: number of input points must'
print,'% be larger than or equal to 4!!!'
print,'% Exiting.....'
radius=-1.
goto,endoffcccode001s
endif
nicc=npoints-2
;
chrdcent_vec=dblarr(2,nicc) ; vector of coordinates of
; the centers of the cir-
; cle chords
dirvec_perpchords=dblarr(2,nicc) ; direction vectors per-
; pendicular to the chords
;
for icc=0L,nicc-1L do begin
chrdcent_vec[0,icc]=(points[0,icc]+points[0,icc+2L])*0.5d0
chrdcent_vec[1,icc]=(points[1,icc]+points[1,icc+2L])*0.5d0
dirvec_perpchords[0,icc]=points[1,icc+2L]-points[1,icc]
dirvec_perpchords[1,icc]=points[0,icc]-points[0,icc+2L]
endfor
;
nicc1=nicc-1L
center_estvec=dblarr(2,nicc1)
for icc=0L,nicc1-1L do begin
m1=dirvec_perpchords[0,icc+1L]*(chrdcent_vec[1,icc]-chrdcent_vec[1,icc+1L])
m2=dirvec_perpchords[1,icc+1L]*(chrdcent_vec[0,icc]-chrdcent_vec[0,icc+1L])
m3=dirvec_perpchords[1,icc+1L]*dirvec_perpchords[0,icc]-dirvec_perpchords[1,icc]*dirvec_perpchords[0,icc+1L]
kk=(m1-m2)/m3
xx=dirvec_perpchords[0,icc]*kk+chrdcent_vec[0,icc]
yy=dirvec_perpchords[1,icc]*kk+chrdcent_vec[1,icc]
center_estvec[0,icc]=xx
center_estvec[1,icc]=yy
endfor
if(nicc1 gt 1) then center=total(center_estvec,2)/double(nicc1) else center=center_estvec
radestvec=dblarr(npoints)
for ip=0L,npoints-1L do begin
val=sqrt((points[0,ip]-center[0])^2+(points[1,ip]-center[1])^2)
radestvec[ip]=val
endfor
radius=total(radestvec)/double(npoints)
;
endoffcccode001s:
end
;
;
;
function monnum2str,month
monstr=' '
month=fix(month)
case month of
1 : monstr='Jan'
2 : monstr='Feb'
3 : monstr='Mar'
4 : monstr='Apr'
5 : monstr='May'
6 : monstr='Jun'
7 : monstr='Jul'
8 : monstr='Aug'
9 : monstr='Sep'
10 : monstr='Oct'
11 : monstr='Nov'
12 : monstr='Dec'
ELSE : monstr=' '
endcase
;
return,monstr
end
;
;
;
pro maxabs_minemis,spectrum,cont_areas,line_area,ratio_maxabs_minemis,tolerance_perc=tolperc
;
if(n_elements(tolperc) eq 0) then tolperc=0.d0 else tolperc=tolperc/100.d0
;
n_pairs=long(float(n_elements(cont_areas))/2.)
xdim=n_elements(spectrum[*,0])
ydim=n_elements(spectrum[0,*])
;
maxabsvec=dblarr(ydim)
minemisvec=dblarr(ydim)
i_emis=0L
i_abs=0L
;
for iy=0L,ydim-1L do begin
avgcval=0.d0
npoints=0L
prof=reform(spectrum[*,iy])
for icnt=0L,n_pairs-1L do begin
cc=2L*icnt
avgcval=avgcval+total(prof[cont_areas[cc]:cont_areas[cc+1L]])
npoints=npoints+n_elements(prof[cont_areas[cc]:cont_areas[cc+1L]])
endfor
avgcval=avgcval/double(npoints)
avglval=total(prof[line_area[0]:line_area[1]])/double(n_elements(prof[line_area[0]:line_area[1]]))
tolint=[avgcval*(1.d0-tolperc),avgcval*(1.d0+tolperc)]
if((avglval lt tolint[0]) or (avglval gt tolint[1])) then begin
if(avglval gt avgcval) then begin ; emission profiles
minemisvec[i_emis]=avgcval
i_emis=i_emis+1L
endif
if(avglval lt avgcval) then begin ; absorption profile
maxabsvec[i_abs]=avgcval
i_abs=i_abs+1
endif
endif
endfor
minemis=0.0d0
maxabs=0.0d0
if(i_emis gt 0) then begin
minemisvec=minemisvec[0:i_emis-1]
minemis=min(minemisvec)
endif
if(i_abs gt 0) then begin
maxabsvec=maxabsvec[0:i_abs-1]
maxabs=max(maxabsvec)
endif
;
;
zerswsw=0
if(minemis le 0.) then zerswsw=zerswsw+1
if(maxabs le 0.) then zerswsw=zerswsw+1
if(zerswsw eq 0) then ratio_maxabs_minemis=minemis/maxabs
;
end
;
;
;
function replace_undefbyunity,mat,valr
; replaces undefined numbers such as NaN and Inf, in matrix mat
; by a value valr. The function can handle scalars, vectors or 2D matrices
;
sss=size(mat)
if(sss[0] eq 0) then begin
t_sw=finite(mat)
if(t_sw eq 0) then result=valr else result=mat
endif
if(sss[0] eq 1) then begin
dim=sss[1]
result=dblarr(dim)
for i=0L,dim-1L do begin
val=mat[i]
t_sw=finite(val)
if(t_sw eq 0) then result[i]=valr else result[i]=val
endfor
endif
if(sss[0] eq 2) then begin
xdim=sss[1]
ydim=sss[2]
result=dblarr(xdim,ydim)
for ix=0L,xdim-1L do begin
for iy=0L,ydim-1L do begin
val=mat[ix,iy]
t_sw=finite(val)
if(t_sw eq 0) then result[ix,iy]=valr else result[ix,iy]=val
endfor
endfor
endif
;
return,result
end
;
;
;
;
;
;
; +++++++++++++++++++++++++++++++++++++++++
; + +
; + MAIN CODE +
; + +
; +++++++++++++++++++++++++++++++++++++++++
;
;
common setup_cb,idl_4degpolyfitting_sw,polyfitting_path
make_setup
if(!version.os_family eq 'unix') then begin
spawn,'pwd',flnmpath
dir_separator='/'
windev='X'
device,retain=2,dec=0
endif
if(!version.os_family eq 'Windows') then begin
spawn,'echo %CD%',flnmpath
dir_separator='\'
windev='win'
endif
;
for itmp=0,4 do print,' '
print,'****************************************************'
print,'* Warning: *'
print,'* This version of the mfsff program must be run *'
print,'* in SolarSoft environment!!!!!!! *'
print,'**************************************************** '
print,' '
print,'Are you running sswidl? (Y/n)'
answ_sswqstn=' '
read,answ_sswqstn
answ_sswqstn=strupcase(answ_sswqstn)
if(answ_sswqstn eq 'N') then begin
print,'Exit IDL and run sswidl'
goto,eocmffffs
endif
for itmp=0,10 do print,' '
;
flnmpath=flnmpath[0]
set_plot,windev
device,retain=2,dec=0
loadct,0
;
mod_proc=0
print,' '
print,' '
print,' +++++++++++++++++++++++++++++++++++++++++'
print,' + +'
print,' + flat-fielding of spectra from MFS +'
print,' + +'
print,' +++++++++++++++++++++++++++++++++++++++++'
answ1=' '
print,' '
print,' MENU:'
print,'==================================================='
print,' 1 .... construct&use flatfield'
print,' 2 .... construct flatfield (default)'
print,' 3 .... use flatfield'
print,' X .... exit'
print,' '
read,answ1
answ1=strupcase(answ1)
case answ1 of
'1' : mod_proc=1
'2' : mod_proc=2
'3' : mod_proc=3
'X' : goto,eocmffffs
ELSE: mod_proc=2
endcase
;
if((mod_proc eq 1) or (mod_proc eq 2)) then begin
print,''
print,'==================================================='
print,' DARK FRAMES'
print,'==================================================='
print,''
answ2=' '
print,'Take all files with ending ''D.fts'' from a specific'
print,' directory? (Y/n)'
read,answ2
answ2=strupcase(answ2)
if(answ2 ne 'N') then answ2='Y'
if(answ2 eq 'Y') then begin
fls_dcd=dialog_pickfile(path='c:\MFS_zpracovani',/directory)
fls_dc=file_search(fls_dcd,'*D.fts',count=n_fls_dc)
endif else begin
fls_dc=dialog_pickfile(path='c:\MFS_zpracovani',/multiple_files,filter=['*D.fts;*d.fts','*.fts','*.*'])
n_fls_dc=n_elements(fls_dc)
endelse
print,' '
print,'% ',n_fls_dc,' files chosen'
fits_read,fls_dc[0],tmpmat1,tmphdr1
xdim=n_elements(tmpmat1[*,0])
ydim=n_elements(tmpmat1[0,*])
dcmat=dblarr(xdim,ydim)*0.d0
for in=0L,n_fls_dc-1L do begin
fits_read,fls_dc[in],tmpmat1,tmphdr1
print,'% reading file ',file_basename(fls_dc[in])
tmpmat1=double(tmpmat1)
dcmat=dcmat+tmpmat1
endfor
dcmat=dcmat/double(n_fls_dc)
print,' '
print,'==================================================='
print,' FLAT FRAMES'
print,'==================================================='
print,' '
print,''
answ3=' '
print,'Take all files with ending ''F.fts'' from a specific'
print,' directory? (Y/n)'
read,answ3
answ3=strupcase(answ3)
if(answ3 ne 'N') then answ3='Y'
if(answ3 eq 'Y') then begin
fls_ffd=dialog_pickfile(path='c:\MFS_zpracovani',/directory)
fls_ff=file_search(fls_ffd,'*F.fts',count=n_fls_ff)
endif else begin
fls_ff=dialog_pickfile(path='c:\MFS_zpracovani',/multiple_files,filter=['*F.fts;*f.fts','*.fts','*.*'])
n_fls_ff=n_elements(fls_ff)
endelse
print,' '
print,'% ',n_fls_ff,' files chosen'
ffmat=dblarr(xdim,ydim)*0.d0
for in=0L,n_fls_ff-1L do begin
print,'% reading file ',file_basename(fls_ff[in])
fits_read,fls_ff[in],tmpmat1,tmphdr1
tmpmat1=double(tmpmat1)
ffmat=ffmat+(tmpmat1-dcmat)
endfor
ffmat=ffmat/double(n_fls_ff)
window,0,xs=xdim,ys=ydim,title='mean FLAT FRAME'
showimg,data=ffmat
xloadct,/block, UPDATECALLBACK='showimg',UPDATECBDATA=ffmat
tvlct,rcol1,gcol1,bcol1,/get
rcol1[254]=0
gcol1[254]=200
bcol1[254]=100
rcol1[255]=240
gcol1[255]=80
bcol1[255]=190
; horizontal wires
print,' '
print,'Click on the two horizontal spectrograph wires'
print,'press enter'
rrr=get_kbrd()
window,0,xs=xdim,ys=ydim,title='mean FLAT FRAME'
showimg,data=ffmat
print,'mark one of the wires'
cursor,xxx,yyy,/normal,/wait
yyy=((yyy>0.)<1.)
yyy_prev=yyy
ypw1=fix(yyy*double(ydim-1L))
print,ypw1,' pxl'
print,'mark other wire'
lggg_thwsp01z:
cursor,xxx,yyy,/normal,/wait
yyy=((yyy>0.)<1.)
if(yyy eq yyy_prev) then goto,lggg_thwsp01z
ypw2=fix(yyy*double(ydim-1L))
print,ypw2,' pxl'
print,'press enter'
rrr=get_kbrd()
ypos_wires=[min([ypw1,ypw2]),max([ypw1,ypw2])]
xpositions=findgen(xdim)
ypositions=findgen(ydim)
;
;
ffmatb=bytscl(ffmat,top=253)
n_tplp=0L
tanhorizincl=0.d0
print,' '
print,' '
print,'=============================-========================'
print,' CORRECTION OF HORIZONTAL INCLINATION OF THE SPECTRUM'
print,'============================-========================='
print,' '
pnts_tmp1=lonarr(2,10000)
tvlct,rcol1,gcol1,bcol1
nxtstrpl01a:
window,0,xs=xdim,ys=ydim,title='mean FLAT FRAME'
tv,ffmatb
contour,ffmatb,xpositions,ypositions,position=[0,0,xdim-1,ydim-1],$
xstyle=1,ystyle=1,/noerase,/nodata,xticklen=-1,yticklen=-1
if(n_tplp gt 0) then begin
for i_tp=0L,n_tplp-1L do begin
plots,pnts_tmp1[0,i_tp],pnts_tmp1[1,i_tp],psym=1,color=254,$
symsize=2.5,thick=2.5,/data
endfor
endif
print,'Mark two points on some horizontal stripes'
print,'The best is to chose one point on the left edge'
print,' and another on the right edge of the spectrum'
print,'!!!! but not on the artificial spectrograph wires'
print,'the first point'
print,' '
print,'press enter'
rrr=get_kbrd()
window,0,xs=xdim,ys=ydim,title='mean FLAT FRAME'
tv,ffmatb
contour,ffmatb,xpositions,ypositions,position=[0,0,xdim-1,ydim-1],$
xstyle=1,ystyle=1,/noerase,/nodata,xticklen=-1,yticklen=-1
if(n_tplp gt 0) then begin
for i_tp=0L,n_tplp-1L do begin
plots,pnts_tmp1[0,i_tp],pnts_tmp1[1,i_tp],psym=1,color=254,$
symsize=2.5,thick=2.5,/data
endfor
endif
print,' '
print,'click left button to mark the first point'
print,'press right mouse button to quit marking'
cursor,xxx1,yyy1,/normal,/wait
if(!mouse.button eq 4) then goto,cntaftrpl
print,fix(xxx1*float(xdim-1L)),' pxl, ',fix(yyy1*float(ydim-1L)),' pxl'
plots,[0,xdim-1L],[fix(yyy1*float(ydim-1L)),fix(yyy1*float(ydim-1L))],linestyle=2,$
color=255,/data,thick=1.5
pnts_tmp1[0,n_tplp]=fix(xxx1*float(xdim-1L))
pnts_tmp1[1,n_tplp]=fix(yyy1*float(ydim-1L))
print,' '
print,' '
print,'click left button to mark the second point'
window,0,xs=xdim,ys=ydim,title='mean FLAT FRAME'
contour,ffmatb,xpositions,ypositions,position=[0,0,xdim-1,ydim-1],$
xstyle=1,ystyle=1,/noerase,/nodata,xticklen=-1,yticklen=-1
tv,ffmatb
if(n_tplp gt 0) then begin
for i_tp=0L,n_tplp-1L do begin
plots,pnts_tmp1[0,i_tp],pnts_tmp1[1,i_tp],psym=1,color=254,$
symsize=2.5,thick=2.5,/data
endfor
endif
plots,[0,xdim-1L],[fix(yyy1*float(ydim-1L)),fix(yyy1*float(ydim-1L))],linestyle=2,$
color=255,/data,thick=1.5
wshow,0
lggg_thwsp02z:
cursor,xxx2,yyy2,/normal,/wait
if((xxx2 eq xxx1) and (yyy2 eq yyy1)) then goto,lggg_thwsp02z
print,fix(xxx2*float(xdim-1L)),' pxl, ',fix(yyy2*float(ydim-1L)),' pxl'
print,' '
print,' '
tanhorizincl=tanhorizincl+(yyy1-yyy2)/(xxx1-xxx2)
n_tplp=n_tplp+1L
goto,nxtstrpl01a
cntaftrpl:
tanhorizincl=tanhorizincl/double(n_tplp)
print,'calculating the inclination angle'
print,'please wait ...'
print,''
print,'tangens of the incl angle: ',tanhorizincl
mm=[min(ffmat),max(ffmat)]
rmvar,pnts_tmp1
;
; excluding horizontal wires and interpolating the intensities there
halfwidth_wires_pxl=8
ffmat_s1=ffmat*0.0d0
wires_yposarea=[ypos_wires[0]-halfwidth_wires_pxl,ypos_wires[0]+halfwidth_wires_pxl,$
ypos_wires[1]-halfwidth_wires_pxl,ypos_wires[1]+halfwidth_wires_pxl]
outside_wires=fltarr(ydim)
iosw=0L
for iy=0L,ydim-1L do begin
sw_noaddprof=0
if((iy ge wires_yposarea[0]) and (iy le wires_yposarea[1])) then sw_noaddprof=sw_noaddprof+1
if((iy ge wires_yposarea[2]) and (iy le wires_yposarea[3])) then sw_noaddprof=sw_noaddprof+1
if(sw_noaddprof eq 0) then begin
outside_wires[iosw]=float(iy)
iosw=iosw+1L
endif
endfor
outside_wires=outside_wires[0:iosw-2L]
for ix=0,xdim-1L do begin
tmpvec1=reform(ffmat[ix,outside_wires])
lin_interp,outside_wires,tmpvec1,ypositions,tmpvec2
ffmat_s1[ix,*]=tmpvec2
endfor
indcsorng=where(ffmat_s1 lt mm[0])
if(indcsorng[0] ne (-1)) then ffmat_s1[indcsorng]=mm[1]
ffmat_s1=(ffmat_s1 print,'% correcting the mean FF for inclination'
print,' '
print,'computations run'
print,'please, wait ...'
print,' '
corr_incl,ffmat_s1,ffmat_s2,tanhorizincl
print,'Done. You can continue ... '
print,' '
;
loadct,0
window,0,xs=xdim,ys=ydim,title='mean FLAT FRAME corrected for inclination'
tvscl,ffmat_s2
;
print,' '
print,' '
print,'=============================-========================'
print,' CORRECTION CURVATURE OF THE SPECTRUM '
print,'============================-========================='
print,' '

print,'Choose area where sp. line used for curvature estimation,'
print,' occurs'
xcurvature=fltarr(2,ydim)*0.
print,' '
print,'press enter'
rrr=get_kbrd()
window,0,xs=xdim,ys=ydim,title='mean FLAT FRAME corrected for inclination'
tvscl,ffmat_s2
print,'mark left lower border'
xxx=0. & yyy=0.
cursor,xxx,yyy,/normal,/wait
xxx=((xxx>0.)<1.)
yyy=((yyy>0.)<1.)
xxx_prev=xxx & yyy_prev=yyy
xbxc1=fix(xxx*float(xdim-1L))
ybxc1=fix(yyy*float(ydim-1L))
print,xbxc1,' pxl, ',ybxc1,' pxl'
print,'mark right upper border'
xxx=0. & yyy=0.
lggg_thwsp03z:
cursor,xxx,yyy,/normal,/wait
if((xxx eq xxx_prev) and (yyy eq yyy_prev)) then goto,lggg_thwsp03z
xxx=((xxx>0.)<1.)
yyy=((yyy>0.)<1.)
xbxc2=fix(xxx*float(xdim-1L))
ybxc2=fix(yyy*float(ydim-1L))
print,xbxc2,' pxl, ',ybxc2,' pxl'
print,' '
print,' '
xbxc=[min([xbxc1,xbxc2]),max([xbxc1,xbxc2])]
ybxc=[min([ybxc1,ybxc2]),max([ybxc1,ybxc2])]
iccv=0L
for iy=0,ydim-1L do begin
if((iy ge ybxc[0]) and (iy le ybxc[1])) then begin
prof_iy=reform(ffmat_s2[xbxc[0]:xbxc[1],iy])
rrr=min(prof_iy,ixrelmin)
xcurvature[0,iccv]=float(iy)
xcurvature[1,iccv]=ixrelmin+xbxc[0]
iccv=iccv+1L
endif
endfor
niccv=iccv-1L
xcurvature=xcurvature[*,0:niccv-1L]
xcurvature[1,*]=xcurvature[1,*]-min(xcurvature[1,*])
xcurvature_sm=xcurvature
window,2
plot,xcurvature[0,*],xcurvature[1,*],xtitle='position along the slit [pxl]',$
psym=1,ytitle='X-shift [pxl]'
crvpsmf=1
shcrvpagl:
read,crvpsmf,prompt='smoothing factor (give 1 for no smoothing): '
xcurvature_sm[1,*]=smooth(xcurvature[1,*],[1,crvpsmf],/edge_truncate)
window,2
plot,xcurvature_sm[0,*],xcurvature_sm[1,*],xtitle='position along the slit [pxl]',$
psym=1,ytitle='X-shift [pxl]'
answcrvp_sw=' '
print,' '
print,'are you satisfied? (Y/n)'
read,answcrvp_sw
answcrvp_sw=strupcase(answcrvp_sw)
if(answcrvp_sw eq 'N') then goto,shcrvpagl
if(idl_4degpolyfitting_sw gt 0) then begin
merrz=dblarr(niccv)+0.5d0
polyfit4deg,xcurvature_sm[0,*],xcurvature_sm[1,*],polyparams,$
meas_err=merrz,winn=3
wdelete,3
endif else begin
openw,1,'ddd.dat'
for idt=0L,niccv-1L do printf,1,xcurvature_sm[0,idt],xcurvature_sm[1,idt]
close,1
spawn,polyfitting_path
polyparams=dblarr(100)
i_polyparams=-1L
val1=0.d0
openr,1,'fitparams.dat'
rdopraggl01:
i_polyparams=i_polyparams+1L
readf,1,val1
polyparams[i_polyparams]=val1
if(eof(1) ne 1) then goto,rdopraggl01
close,1
polyparams=polyparams[0:i_polyparams]
n_polyparams=n_elements(polyparams)
spawn,'rm -f fitparams.dat ddd.dat'
print,''
endelse
print,''
print,'% correcting the mean FF for curvature'
print,' '
print,'computations run'
print,'please, wait ...'
print,' '
corr_curv,ffmat_s2,ffmat_s3,polyparams
print,'Done. You can continue ... '
print,' '
print,'press enter'
window,0,xs=xdim,ys=ydim,title='mean FLAT FRAME corrected for both inclination & curvature'
tvscl,ffmat_s3
rrr=get_kbrd()
;
loadct,0
print,' '
wdelete,2
print,' '
print,' '
print,' '
print,'=============================-========================'
print,' SOFT-FLAT MATRIX '
print,'============================-========================='
print,' '
print,'% smoothing the corrected mean FF'
print,' '
print,'computations run'
print,'please, wait ...'
print,' '
tmpmat1=smooth(ffmat_s3,[9,9],/edge_truncate)
ffmat_ms1=smooth(tmpmat1,[1,9],/edge_truncate)
;
ffmat_ms2=dblarr(xdim,ydim)
xposvec_intbox=ysection(min(xpositions),max(xpositions),3.)
for iy=0L,ydim-1L do begin
tmpvec1=reform(ffmat_ms1[xposvec_intbox,iy])
tmpvec2=interpol(tmpvec1,xposvec_intbox,xpositions,/LSQUADRATIC)
ffmat_ms2[*,iy]=tmpvec2
endfor
;
yposvec_intbox=ysection(min(ypositions),max(ypositions),9.)
ffmat_ms3=dblarr(xdim,ydim)
for ix=0L,xdim-1L do begin
tmpvec1=reform(ffmat_ms2[ix,yposvec_intbox])
lin_interp,yposvec_intbox,tmpvec1,ypositions,tmpvec2
ffmat_ms3[ix,*]=tmpvec2
endfor
print,'Done.'
print,' '
;
; removing spectral lines from the soft flat
print,'% removing spectral lines from the soft flat'
avgprof_ms=dblarr(xdim)*0.0d0
for iy=0L,ydim-1L do avgprof_ms=avgprof_ms+reform(ffmat_ms3[*,iy])
avgprof_ms=avgprof_ms/double(ydim)
indcs_zer=where(avgprof_ms eq 0)
if(indcs_zer[0] ne (-1)) then avgprof_ms[indcs_zer]=1.d-10
ffmat_ms=dblarr(xdim,ydim)
for iy=0L,ydim-1L do begin
ffmat_ms[*,iy]=reform(ffmat_ms3[*,iy])/avgprof_ms
endfor
indcs_zer=(ffmat_ms3 eq 0.)
if(indcs_zer[0] ne (-1)) then ffmat_ms3[indcs_zer]=1.
window,0,xs=xdim,ys=ydim,title='soft flat matrix'
tvscl,hist_filter(ffmat_ms,percentage=1.)
;
;
loadct,0
print,' '
print,' '
print,' '
print,'=============================-========================'
print,' HARD-FLAT MATRIX '
print,'============================-========================='
print,' '
ffmat_mh1=ffmat_s3/ffmat_ms3
print,'% removing any traces of spectral lines from the hard flat'
print,'% if there are still any'
ffmat_mh=dblarr(xdim,ydim)
avgprof_mh=dblarr(xdim)*0.0d0
for iy=0L,ydim-1L do avgprof_mh=avgprof_mh+reform(ffmat_mh1[*,iy])
avgprof_mh=avgprof_mh/double(ydim)
indcs_zer=where(avgprof_mh le 0.)
if(indcs_zer[0] ne (-1)) then avgprof_mh[indcs_zer]=1.
for iy=0L,ydim-1L do begin
ffmat_mh[*,iy]=reform(ffmat_mh1[*,iy])/avgprof_mh
endfor
loadct,0
window,1,xs=xdim,ys=ydim,title='hard flat matrix'
tvscl,dtfilter(ffmat_mh,threshold=2.)
print,' '
print,'press enter'
rrr=get_kbrd()
print,' '
print,' '
print,' '
print,'=============================-========================'
print,' SLIT-FLAT MATRIX '
print,'============================-========================='
print,' '
ffmat_slit=dblarr(xdim,ydim)
tmpvec1=total(ffmat_mh,1)/double(xdim)
indcs_zer=where(tmpvec1 le 0.)
if(indcs_zer[0] ne (-1)) then tmpvec1[indcs_zer]=1.
for ix=0l,xdim-1L do ffmat_slit[ix,*]=tmpvec1
ffmat_slit=replace_undefbyunity(ffmat_slit,1.)
ffmat_slit=dtfilter(ffmat_slit,threshold=10.)
ffmat_mhp=ffmat_mh/ffmat_slit
ffmat_mhp=dtfilter(ffmat_mhp,threshold=10.)
loadct,0
window,2,xs=xdim,ys=ydim,title='slit-flat matrix (drifting)'
tvscl,dtfilter(ffmat_slit,threshold=2.)
wset,1
tvscl,dtfilter(ffmat_mhp,threshold=2.)
print,' '
print,'press enter'
rrr=get_kbrd()
print,' '
print,' '
print,' '
print,'=============================-========================'
print,' RESULTING FLATFIELD MATRICES '
print,'============================-========================='
print,' '
wdelete,1
loadct,0
ffmat_static=ffmat_mhp*ffmat_ms
ffmat_static=replace_undefbyunity(ffmat_static,1.)
indcs_zer=where(ffmat_static le 0.)
if(indcs_zer[0] ne (-1)) then ffmat_static[indcs_zer]=1.
window,0,xs=xdim,ys=ydim,title='static-flat matrix'
tvscl,hist_filter(ffmat_static,percentage=1.)
;
outffdir=flnmpath
print,'Save flat-field matrices into idl-save file? (y/N)'
answ_fffs1=' '
read,answ_fffs1
answ_fffs1=strupcase(answ_fffs1)
if(answ_fffs1 eq 'Y') then begin
answ_fffs2=''
print,' '
print,'The file is to be saved into directory:'
print,outffdir
print,'Do you wanna change this path? (y/N)'
read,answ_fffs2
answ_fffs2=strupcase(answ_fffs2)
if(answ_fffs2 eq 'Y') then begin
outffdir=dialog_pickfile(path='c:\MFS_zpracovani',/directory,title='choose the directory')
endif
fffilename=''
read,fffilename,prompt='filename without extension: '
ll1z=strlen(outffdir)
ll2z=strmid(outffdir,ll1z-1L,1)
if(ll2z ne dir_separator) then outffdir=outffdir+dir_separator
fffilename=outffdir+fffilename+'.idl'
save,ffmat_static,ffmat_slit,tanhorizincl,ypos_wires,$
polyparams,dcmat,file=fffilename
print,' '
print,'Flat-field matrices saved into:'
print,fffilename
print,' '
endif
wdelete,0
wdelete,2
;
endif
;
if(mod_proc eq 3) then begin
flinff=dialog_pickfile(path='c:\MFS_zpracovani',title='idl-save file with FF matrices')
print,' '
print,'reading file:'
print,flinff
restore,flinff,/v
xdim=n_elements(ffmat_static[*,0])
ydim=n_elements(ffmat_static[0,*])
ypositions=findgen(ydim)
xpositions=findgen(xdim)
print,' '
endif
;
;
if((mod_proc eq 1) or (mod_proc eq 3)) then begin
print,' '
print,' '
print,' '
print,'=============================-========================'
print,' FLAT-FIELDING '
print,'============================-========================='
print,' '
flinsp=dialog_pickfile(path='c:\MFS_zpracovani',title='fts file with spectrum')
print,' '
print,'reading file:'
print,flinsp
fits_read,flinsp,spmat,sphdr
sep_pathflnm,flinsp,path_sp,flnm_sp,flext_sp
spmat=double(spmat)
xpp=findgen(xdim)
window,0,xs=xdim,ys=ydim,title='raw spectra'
showimg,data=spmat
xloadct,/block, UPDATECALLBACK='showimg',UPDATECBDATA=spmat
; horizontal wires
print,' '
print,'Click on the two horizontal spectrograph wires'
print,'press enter'
rrr=get_kbrd()
window,0,xs=xdim,ys=ydim,title='raw spectra'
showimg,data=spmat
print,'mark one of the wires'
cursor,xxx,yyy,/normal,/wait
yyy=((yyy>0.)<1.)
yyy_prev=yyy
ypw1=fix(yyy*double(ydim-1L))
print,ypw1,' pxl'
print,' '
print,'mark other wire'
lggg5_z01abcd:
cursor,xxx,yyy,/normal,/wait
yyy=((yyy>0.)<1.)
if(yyy eq yyy_prev) then goto,lggg5_z01abcd
ypw2=fix(yyy*double(ydim-1L))
print,ypw2,' pxl'
print,'press enter'
rrr=get_kbrd()
ypos_wires_sp=[min([ypw1,ypw2]),max([ypw1,ypw2])]
print,'% subtracting dark-frame'
spmat1=spmat-dcmat
print,'% correcting the spectrum for inclination'
print,' '
print,'computations run'
print,'please, wait ...'
print,' '
corr_incl,spmat1,spmat2,tanhorizincl
print,'Done.'
print,' '
print,'% correcting the spectrum for curvature'
print,' '
print,'computations run'
print,'please, wait ...'
print,' '
corr_curv,spmat2,spmat3,polyparams
print,'Done. '
spmat3=replace_undefbyunity(spmat3,1.)
spmat3=(spmat3>0.)
loadct,0
window,0,xs=xdim,ys=ydim,title='corrected raw spectra'
tvscl,spmat3
print,' '
print,' '
print,'In the raw corrected spectra chose area at the disc'
print,' to be used for estimation of vertical shift of the slit-flat frame'
print,' '
print,'click into image to fix the first border of the area'
cursor,xxx,yyy,/wait,/normal
yyy=((yyy>0.)<1.)
yyy_prev=yyy
y1_avgp=fix(yyy*float(ydim-1L))
print,y1_avgp,' pxl'
print,' '
print,'click into image to fix other border of the area'
lggg5_z02abcd:
cursor,xxx,yyy,/wait,/normal
yyy=((yyy>0.)<1.)
if(yyy eq yyy_prev) then goto,lggg5_z02abcd
y2_avgp=fix(yyy*float(ydim-1L))
print,y2_avgp,' pxl'
y1_avgp1=min([y1_avgp,y2_avgp])
y2_avgp2=max([y1_avgp,y2_avgp])
sp_ffldbystatff=spmat3/ffmat_static
estim_ffshift,sp_ffldbystatff[*,y1_avgp1:y2_avgp2],ffmat_slit[*,y1_avgp1:y2_avgp2],ffydrift,$
shift_range=[-30.,30.],shift_step=1.,/auto,pltwindow=5
sydrft=string(ffydrift,format='(f7.2)')
print,' '
print,' '
print,'Y drift of '+sydrft+' pxl of the slit flat matrix was found automatically.'
print,'Are you satisfied? (y/N)'
ydrft_sw0=' '
read, ydrft_sw0
ydrft_sw0=strupcase(ydrft_sw0)
if(ydrft_sw0 ne 'Y') then $
estim_ffshift,sp_ffldbystatff[*,y1_avgp1:y2_avgp2],ffmat_slit[*,y1_avgp1:y2_avgp2],ffydrift,$
pltwindow=5
rmvar,sp_ffldbystatff
print,' '
print,'Apply the drift? (y/N)'
ydrft_sw=' '
read, ydrft_sw
wdelete,5
print,' '
ydrft_sw=strupcase(ydrft_sw)
if(ydrft_sw eq 'Y') then begin
ffmat_slit1=dblarr(xdim,ydim)*0.d0
yposshftd=ypositions-ffydrift
print,'% shifting the slit-flat matrix'
print,' '
print,'computations run'
print,'please, wait ...'
print,' '
for ix=0L,xdim-1L do begin
tmpvec1=ffmat_slit[ix,*]
lin_interp,yposshftd,tmpvec1,ypositions,tmpvec2,/outside_const
ffmat_slit1[ix,*]=tmpvec2
endfor
print,' '
print,'Done.'
endif else begin
ffmat_slit1=ffmat_slit
endelse
fltfield=ffmat_static*ffmat_slit1
min_fff=min(fltfield)
if(min_fff lt 0.) then fltfield=fltfield-min_fff ; save FF from negative and zero values
indcs_zer=where(fltfield eq 0.)
if(indcs_zer[0] ne (-1)) then fltfield[indcs_zer]=1.
avgffval=total(fltfield)/double(xdim*ydim)
fltfield=((fltfield/avgffval>0.1)<10.) ; final flatfield
ii_zero=where(fltfield eq 0.)
if(ii_zero[0] ne (-1)) then fltfield[ii_zero]=1. ; replace zeroes by 1.
spflfld=spmat3/fltfield ; flatfielded spectrum
window,1,xs=xdim,ys=ydim,title='flatfielded raw spectra'
tvscl,spflfld
window,2,xs=xdim,ys=ydim,title='flatfield'
tvscl,hist_filter(fltfield)
print,' '
print,' '
print,'press enter'
rrr=get_kbrd()
wdelete,0
wdelete,2
wset,1
avg_obs_prof_mat=spflfld[*,y1_avgp1:y2_avgp2]
n_avg_obs_prof_mat=n_elements(avg_obs_prof_mat[0,*])
avg_obs_prof=total(avg_obs_prof_mat,2)/float(n_avg_obs_prof_mat)
rmvar,avg_obs_prof_mat
rmvar,n_avg_obs_prof_mat
;
restore,'c:\MFS_zpracovani\software\mHalp_samplespectra.idl'

window,0,title='sample avg disc spectrum [ window No.1 ]'
plot,SHAAVGPROF[0,*],SHAAVGPROF[1,*],xtitle='!7k!3 [!z(c5)]',ytitle='arbitrary units'
window,2,title='avg observed spectrum [ window No.2 ]'
plot,xpositions,avg_obs_prof/max(avg_obs_prof),xtitle='!3wavelength pixel',ytitle='arbitrary units'
print,' '
print,'Compare profiles in the windows No.1 and No.2'
print,'Does it seem that the observed spectrum is reversed in wavelength? (y/N)'
answ_rvsdwvl_sw1=' '
read,answ_rvsdwvl_sw1
print,' '
answ_rvsdwvl_sw1=strupcase(answ_rvsdwvl_sw1)
if(answ_rvsdwvl_sw1 eq 'Y') then begin
spflfld1=reverse(spflfld,1)
window,1,xs=xdim,ys=ydim,title='flatfielded raw spectra'
tvscl,spflfld1
window,0,title='sample avg disc spectrum'
plot,SHAAVGPROF[0,*],SHAAVGPROF[1,*],xtitle='!7k!3 [!z(c5)]',ytitle='arbitrary units'
window,2,title='avg observed spectrum reversed in wvl'
plot,xpositions,reverse(avg_obs_prof)/max(avg_obs_prof),xtitle='!3wavelength pixel',ytitle='arbitrary units'
print,'Are you satisfied? (Y/n)'
answ_rvsdwvl_sw2=' '
read,answ_rvsdwvl_sw2
answ_rvsdwvl_sw2=strupcase(answ_rvsdwvl_sw2)
if(answ_rvsdwvl_sw2 ne 'N') then begin
spflfld=temporary(spflfld1)
endif else begin
window,2,title='avg observed spectrum reversed in wvl'
plot,xpositions,avg_obs_prof/max(avg_obs_prof),xtitle='!3wavelength pixel',ytitle='arbitrary units'
window,1,xs=xdim,ys=ydim,title='flatfielded raw spectra'
tvscl,spflfld
rmvar,spflfld1
endelse
endif
print,' '
print,' '
sw_wvl=' '
print,'Estimate dispersion and create wvl vector? (Y/n)'
read,sw_wvl
sw_wvl=strupcase(sw_wvl)
if(sw_wvl ne 'N') then sw_wvl='Y'
if (sw_wvl eq 'Y') then begin
rcol=[0.,1.,1.,0.]
gcol=[0.,1.,0.,1.]
bcol=[0.,1.,0.,0.]
tvlct,rcol*255,gcol*255,bcol*255
avg_obs_prof_mat=spflfld[*,y1_avgp1:y2_avgp2]
n_avg_obs_prof_mat=n_elements(avg_obs_prof_mat[0,*])
avg_obs_prof=total(avg_obs_prof_mat,2)/float(n_avg_obs_prof_mat)
rmvar,avg_obs_prof_mat
rmvar,n_avg_obs_prof_mat
window,0,title='sample avg disc spectrum'
plot,SHAAVGPROF[0,*],SHAAVGPROF[1,*],xtitle='!7k!3 [!z(c5)]',ytitle='arbitrary units'
window,2,title='avg observed spectrum'
avg_obs_prof_norm=avg_obs_prof/max(avg_obs_prof)
plot,xpositions,avg_obs_prof_norm,xtitle='!3wavelength pixel',ytitle='arbitrary units'
print,' '
print,'mark at least two points in avg observed spectrum'
print,' press right mouse-button to quit marking'
print,' '
wvl_points=fltarr(1000)
pxl_points=fltarr(1000)
iwp=-1
xxx_orev=-1
nxtwptl:
iwp=iwp+1
lggg_wvlest:
cursor,xxx,yyy,/wait,/data
if((!mouse.button eq 4) and (iwp ge 1)) then goto,no_nxtwptl
if(xxx eq xxx_orev) then goto,lggg_wvlest else xxx_orev=xxx
print,iwp+1,' ',xxx,' pxl'
rrr=min(abs(xxx-xpositions),ipos)
oplot,[xxx],[avg_obs_prof_norm[ipos]],color=2,psym=1,symsize=2
xyouts,xxx,avg_obs_prof_norm[ipos]-0.05,strcompress(string(iwp+1),/remove_all)
pxl_points[iwp]=xxx
goto,nxtwptl
no_nxtwptl:
nwp=iwp
print,' '
print,'Press enter'
rrr=get_kbrd()
wshow,2
print,' '
print,'Now, mark the same wavelength points in the sample spectrum'
window,0,title='sample avg disc spectrum'
plot,SHAAVGPROF[0,*],SHAAVGPROF[1,*],xtitle='!7k!3 [!z(c5)]',ytitle='arbitrary units'
xxx_orev=-1
for iwp=0,nwp-1 do begin
lggg_wvlest01:
cursor,xxx,yyy,/wait,/data
if(xxx eq xxx_orev) then goto,lggg_wvlest01 else xxx_orev=xxx
print,iwp+1,' ',xxx,' A'
wvl_points[iwp]=xxx
oplot,[xxx],[yyy],color=2,psym=1,symsize=2
endfor
wdelete,0
wdelete,2
wvl_points=wvl_points[0:nwp-1]
pxl_points=pxl_points[0:nwp-1]
pxl2wvl=linfit(pxl_points,wvl_points,/double)
wvl=pxl2wvl[0]+pxl2wvl[1]*xpositions
print,' '
print,'wvl_o: ',pxl2wvl[0]
print,'dispersion: ',pxl2wvl[1]
print,' '
print,'press enter'
window,0,title='calibration curve for wavelength'
plot,xpositions,wvl,/nodata,ystyle=16,xtitle='!3wvl pixel',ytitle='!3wavelength [!z(c5)]',$
title='calibration curve for wavelength'
oplot,pxl_points,wvl_points,psym=1,color=3
oplot,xpositions,wvl
rrr=get_kbrd()
wdelete,0
endif
loadct,0
print,' '
print,' '
posd1=strpos(flnm_sp,'_')
sd1=strmid(flnm_sp,posd1+1,2)
sd2=strmid(flnm_sp,posd1+3,2)
sd3=strmid(flnm_sp,posd1+5,2)
st1=strmid(flnm_sp,posd1+8,2)
st2=strmid(flnm_sp,posd1+10,2)
st3=strmid(flnm_sp,posd1+12,2)
tst1=strmid(flnm_sp,posd1+14,1)
tst1=strupcase(tst1)
tst1t=where(tst1 eq ['N','F','D','.'])
if(tst1t[0] eq (-1)) then begin
st4=strmid(flnm_sp,posd1+14,2)
endif else begin
st4='00'
endelse
if(fix(sd1) gt 50) then sd1a='19'+sd1 else sd1a='20'+sd1
;
sw_abscal=' '
if(sw_wvl eq 'Y') then begin
print,'Make absolute calibration according to local (close) disc intensity? (Y/n)'
read,sw_abscal
sw_abscal=strupcase(sw_abscal)
endif else begin
sw_abscal='N'
endelse
if(sw_abscal ne 'N') then sw_abscal='Y'
if(sw_abscal eq 'Y') then begin
zdgflt_l001:
window,1,xs=xdim,ys=ydim,title='flatfielded raw spectra'
tvscl,spflfld
answ_advv_sw='N'
print,'please, chose Y borders of the disk in the spectrum image'
print,' '
print,'press enter'
rrr=get_kbrd()
window,1,xs=xdim,ys=ydim,title='flatfielded raw spectra'
tvscl,spflfld
print,' '
print,'click into image to fix the first border'
cursor,xxx,yyy,/wait,/normal
yyy=((yyy>0.)<1.)
yyy_prev=yyy
y1_zfflt=fix(yyy*float(ydim-1L))
print,y1_zfflt,' pxl'
print,' '
print,'click into image to fix other border'
lggg5_z03abcd:
cursor,xxx,yyy,/wait,/normal
yyy=((yyy>0.)<1.)
if(yyy eq yyy_prev) then goto,lggg5_z03abcd
y2_zfflt=fix(yyy*float(ydim-1L))
print,y2_zfflt,' pxl'
y1_spdisk=min([y1_zfflt,y2_zfflt])
y2_spdisk=max([y1_zfflt,y2_zfflt])
spflfld_counts=spflfld
print,' '
print,' '
print,'Pick the slit-jaw complementary to the spectra'
print,' read from the file '+flnm_sp+'.'+flext_sp
fl_sj=dialog_pickfile(path='c:\MFS_zpracovani',title='Pick the slit-jaw fits file')
fits_read,fl_sj,imgsj,hdrsj
sss=size(imgsj)
xdim_sj=sss[1]
ydim_sj=sss[2]
date_sj=strmid(hdrsj[7],11,9)
tmpstr=strmid(hdrsj[8],11)
posap=strpos(tmpstr,'''')
time_sj=strmid(tmpstr,0,posap)
rmvar,tmpstr
rmvar,posap
loadct,0
tvlct,rcol,gcol,bcol,/get
rcol[255]=255
gcol[255]=255
bcol[255]=255
rcol[254]=255
gcol[254]=0
bcol[254]=0
rcol[253]=200
gcol[253]=200
bcol[253]=0
tvlct,rcol,gcol,bcol
imgsjb=bytscl(imgsj,top=252)
estcc_l001:
window,4,xs=xdim_sj,ys=ydim_sj,title='slit-jaw image'
tv,imgsjb
contour,imgsjb,findgen(xdim_sj),findgen(ydim_sj),$
position=[0,0,xdim_sj-1,ydim_sj-1],xstyle=1,ystyle=1,/noerase,/nodata
ypos_wires_sj=dblarr(2,2) ; [x,y] positions of upper and lower hairs in the slit jaw
print,'click on the intersection of upper wire with the slit in the slit jaw'
cursor,xx,yy,/wait,/data
print,xx,' pxl ',yy,' pxl'
ypos_wires_sj[0,1]=xx
ypos_wires_sj[1,1]=yy
xx_prev=xx
yy_prev=yy
print,'click on the intersection of lower wire with the slit in the slit jaw'
lggg_thwsp04z:
cursor,xx,yy,/wait,/data
if((xx eq xx_prev) and (yy eq yy_prev)) then goto,lggg_thwsp04z
print,xx,' pxl ',yy,' pxl'
ypos_wires_sj[0,0]=xx
ypos_wires_sj[1,0]=yy
window,4,xs=xdim_sj,ys=ydim_sj,title='slit-jaw image'
tv,imgsjb
contour,imgsjb,findgen(xdim_sj),findgen(ydim_sj),ticklen=-1,$
position=[0,0,xdim_sj-1,ydim_sj-1],xstyle=1,ystyle=1,/noerase,/nodata
print,' '
print,'mark points (at least 4 but not more than 100) at the limb'
print,'to estimate center and radius of the disc'
print,' '
print,'press enter'
rrr=get_kbrd()
window,4,xs=xdim_sj,ys=ydim_sj,title='slit-jaw image'
tv,imgsjb
contour,imgsjb,findgen(xdim_sj),findgen(ydim_sj),ticklen=-1,$
position=[0,0,xdim_sj-1,ydim_sj-1],xstyle=1,ystyle=1,/noerase,/nodata
limbpoints=fltarr(2,100)
icc=-1L
xx_prev=-1 & yy_prev=-1
chsnxtlmbpt:
print,' '
print,'mark',icc+2,' point on the limb'
if(icc ge 2) then print,'or press right mouse button to quit marking'
lggg_lmbest:
cursor,xx,yy,/wait,/data
if(icc ge 2) then begin
if(!mouse.button eq 4) then goto,nnxtlmbp0z
endif
if((xx eq xx_prev) and (yy eq yy_prev)) then begin
goto,lggg_lmbest
endif else begin
xx_prev=xx
yy_prev=yy
endelse
print,xx,' pxl ',yy,' pxl'
icc=icc+1L
limbpoints[0,icc]=xx
limbpoints[1,icc]=yy
oplot,[xx],[yy],psym=1,symsize=2,color=254
goto,chsnxtlmbpt
nnxtlmbp0z:
limbpoints=limbpoints[*,0:icc]
print,' '
find_circcent,limbpoints,disccenter_pxl,radius=discradius_pxl
npts_limbfit=100L
phivec=findgen(npts_limbfit+1L)/float(npts_limbfit)*2.d0*!dpi
phivec=double(phivec)
limbfit=dblarr(2,npts_limbfit+1L)
limbfit[0,*]=cos(phivec)*double(discradius_pxl)+disccenter_pxl[0]
limbfit[1,*]=sin(phivec)*double(discradius_pxl)+disccenter_pxl[1]
oplot,limbfit[0,*],limbfit[1,*],color=253,thick=2.0
answ_ccok=' '
print,'Are you satisfied? (y/N)'
read,answ_ccok
answ_ccok=strupcase(answ_ccok)
if(answ_ccok ne 'Y') then begin
print,' '
print,'Try different distribution of points at the limb.'
print,' '
goto,estcc_l001
endif
print,' '
print,'disc center in pxl: ',disccenter_pxl
print,'disc radius in pxl: ',discradius_pxl
print,' '
loadct,0
window,1,xs=xdim,ys=ydim,title='flatfielded raw spectra -- real relative intensities'
tvscl,spflfld
print,' '
print,'please, fix place at the disc close observed object, but not too'
print,' close to the limb (mu>0.14) for calibration according a tabulated'
print,' profile'
print,' Profile from this place is to be used for absolute calibration'
print,' '
cursor,xxx,yyy,/wait,/normal
yyy=((yyy>0.)<1.)
ysp_4calib=fix(yyy*float(ydim-1L))
print,ysp_4calib,' pxl in spectrum'
;
; position in the SJ corresponding to ysp_4calib
t_par=(ysp_4calib-ypos_wires_sp[0])/(ypos_wires_sp[1]-ypos_wires_sp[0]) ; parameter of section
; between the wires both
; in spectrum and SJ
xsj_4calib=(ypos_wires_sj[0,1]-ypos_wires_sj[0,0])*t_par+ ypos_wires_sj[0,0]
ysj_4calib=(ypos_wires_sj[1,1]-ypos_wires_sj[1,0])*t_par+ ypos_wires_sj[1,0]
print,' '
print,'corresponding position in slit-jaw in pxl'
print,xsj_4calib,', ',ysj_4calib
print,' '
print,' '
cmpltdtstr=sd1a+sd2+sd3+' '+st1+st2+st3+'.'+st4
rrr=PB0R(cmpltdtstr,/earth,/arcsec)
discradius_arcsec=rrr[2]
print,'disc radius in pxl: ',discradius_pxl
print,'disc radius in arcsec: ',discradius_arcsec
one_pxl_in_arcsec=discradius_arcsec/discradius_pxl
print,'pixel size: ',one_pxl_in_arcsec,' arcsec'
print,' '
xsj_4calib_pxl=xsj_4calib-disccenter_pxl[0]
ysj_4calib_pxl=ysj_4calib-disccenter_pxl[1]
xsj_4calib_arcsec=xsj_4calib_pxl*one_pxl_in_arcsec
ysj_4calib_arcsec=ysj_4calib_pxl*one_pxl_in_arcsec*(-1.d0) ; multiplication by -1 because the MFS observations (both spectra and SJ)
; are upside-down
print,' '
print,'position in ortogonal disc-coordinates in arcsec of the place used for the calibration:'
print,xsj_4calib_arcsec,' ',ysj_4calib_arcsec
print,' '
; heliographic coordinates of the position used for the calibration
bl_4calib=ARCMIN2HEL(xsj_4calib_arcsec/60.,ysj_4calib_arcsec/60.,date=cmpltdtstr)*!dpi/180.d0
bl_0=ARCMIN2HEL(0.,0.,date=cmpltdtstr)*!dpi/180.d0 ; heliographic coordinates of the disc center
muest=sin(bl_4calib[0])*sin(bl_0[0])+cos(bl_4calib[0])*cos(bl_0[0])*cos(bl_4calib[1]-bl_0[1])
print,'mu estimated for the selected position: ',muest
;
wset,1
ypos_sctlofflimb=lonarr(2)
print,' '
print,'fix two off-limb Y-positions outside the prominence for estimation'
print,' of the off-limb scattered light along the slit'
print,'press enter'
rrr=get_kbrd()
window,1,xs=xdim,ys=ydim,title='flatfielded raw spectra -- real relative intensities'
tvscl,spflfld
print,'the first position ...'
cursor,xxx,yyy,/wait,/normal
yyy=((yyy>0.)<1.)
yyy_prev=yyy
y1sctl=fix(yyy*float(ydim-1L))
print,y1sctl,' pxl'
print,' '
print,'other position ... '
lggg_z01ab34:
cursor,xxx,yyy,/wait,/normal
if(yyy eq yyy_prev) then goto,lggg_z01ab34
yyy=((yyy>0.)<1.)
y2sctl=fix(yyy*float(ydim-1L))
print,y2sctl,' pxl'
print,' '
print,' '
ypos_sctlofflimb[0]=min([y1sctl,y2sctl])
ypos_sctlofflimb[1]=max([y1sctl,y2sctl])
nr_aop4c=5
ylsp_4calib=((ysp_4calib-nr_aop4c)>0)
yusp_4calib=((ysp_4calib+nr_aop4c)<(ydim-1L))
avgobsprof4calib=total(spflfld[*,ylsp_4calib:yusp_4calib],2)$
/double(yusp_4calib-ylsp_4calib+1L)
print,' '
print,'press enter'
rrr=get_kbrd()
window,2,title='observed disc profile'
plot,wvl,avgobsprof4calib,xtitle='!7k!3 [!z(c5)]',ytitle='!6I!3[counts]'
print,'click at the profile center'
cursor,xxx,yyy,/wait,/data
print,xxx,' A'
wvlcen_o=xxx
relwvl_o=wvl-wvlcen_o
davidhalpprof,muest,wvl,tabint,relwvl=relwvl_t
lin_interp,relwvl_t,tabint,relwvl_o,tabint_to
indcs_zer=where(avgobsprof4calib le 0)
if(indcs_zer[0] ne (-1)) then avgobsprof4calib[indcs_zer]=1.
calibvec1=tabint_to/avgobsprof4calib
calcoef1=total(calibvec1)/double(n_elements(calibvec1))
print,' '
print,' '
;
; choosing sections along the dispersion of the line and of continum
; on the left and right from the line
rcol2=[0.,1.,1.,0.]
gcol2=[0.,1.,0.,1.]
bcol2=[0.,1.,0.,0.]
tvlct,rcol2*255,gcol2*255,bcol2*255
twocontareas=lonarr(4)
linearea=lonarr(2)
window,2,title='observed disc profile versus spectral pxls'
plot,avgobsprof4calib,xtitle='!7k!3 [pxl]',ytitle='!6I!3[counts]'
mmaop4c=max(avgobsprof4calib)*2.d0
print,'chose continnum area on the left from the spectr. line'
print,'press enter'
r=get_kbrd()
window,2,title='observed disc profile versus spectral pxls'
plot,avgobsprof4calib,xtitle='!7k!3 [pxl]',ytitle='!6I!3[counts]'
print,'left border'
cursor,xxx,yyy,/data,/wait
twocontareas[0]=xxx
xxx_prev=xxx
oplot,[xxx,xxx],[0.,mmaop4c],color=2,line=1
print,xxx,' pxl'
print,' '
print,'right border'
lggg_z01:
cursor,xxx,yyy,/data,/wait
if(xxx eq xxx_prev) then goto,lggg_z01
twocontareas[1]=xxx
oplot,[xxx,xxx],[0.,mmaop4c],color=2,line=1
print,xxx,' pxl'
print,' '
if(twocontareas[0] gt twocontareas[1]) then begin
tmpval=twocontareas[1]
twocontareas[1]=twocontareas[0]
twocontareas[0]=tmpval
endif
print,'chose continnum area on the right from the spectr. line'
print,'press enter'
r=get_kbrd()
window,2,title='observed disc profile versus spectral pxls'
plot,avgobsprof4calib,xtitle='!7k!3 [pxl]',ytitle='!6I!3[counts]'
oplot,[twocontareas[0],twocontareas[0]],[0.,mmaop4c],color=2,line=1
oplot,[twocontareas[1],twocontareas[1]],[0.,mmaop4c],color=2,line=1
print,'left border'
cursor,xxx,yyy,/data,/wait
twocontareas[2]=xxx
xxx_prev=xxx
oplot,[xxx,xxx],[0.,mmaop4c],color=2,line=1
print,xxx,' pxl'
print,' '
print,'right border'
lggg_z02:
cursor,xxx,yyy,/data,/wait
if(xxx eq xxx_prev) then goto,lggg_z02
twocontareas[3]=xxx
oplot,[xxx,xxx],[0.,mmaop4c],color=2,line=1
print,xxx,' pxl'
if(twocontareas[2] gt twocontareas[3]) then begin
tmpval=twocontareas[3]
twocontareas[3]=twocontareas[2]
twocontareas[2]=tmpval
endif
print,' '
print,'chose area of the spectr. line'
print,'press enter'
r=get_kbrd()
window,2,title='observed disc profile versus spectral pxls'
plot,avgobsprof4calib,xtitle='!7k!3 [pxl]',ytitle='!6I!3[counts]'
oplot,[twocontareas[0],twocontareas[0]],[0.,mmaop4c],color=2,line=1
oplot,[twocontareas[1],twocontareas[1]],[0.,mmaop4c],color=2,line=1
oplot,[twocontareas[2],twocontareas[2]],[0.,mmaop4c],color=2,line=1
oplot,[twocontareas[3],twocontareas[3]],[0.,mmaop4c],color=2,line=1
print,'left border'
cursor,xxx,yyy,/data,/wait
linearea[0]=xxx
xxx_prev=xxx
oplot,[xxx,xxx],[0.,mmaop4c],color=2,line=1
print,xxx,' pxl'
print,' '
print,'right border'
lggg_z03:
cursor,xxx,yyy,/data,/wait
if(xxx eq xxx_prev) then goto,lggg_z03
linearea[1]=xxx
oplot,[xxx,xxx],[0.,mmaop4c],color=2,line=1
print,xxx,' pxl'
print,' '
maxabs_minemis,spflfld,twocontareas,linearea,ratio_maxabs_minemis,$
tolerance_perc=8.
print,'a ratio between the continuum intensities at the faintest'
print,' prominence and the brightest disk profiles: ',ratio_maxabs_minemis
print,' '
print,' '
wdelete,2
window,3,title='tabulated profile'
plot,relwvl_o,avgobsprof4calib*calcoef1,xtitle='!7Dk!3 [!z(c5)]',ytitle='!6I!3[erg/cm!u2!n/s/sr/Hz]'
oplot,relwvl_o,tabint_to,color=2
iwp=-1
calibratios=dblarr(4,2000)
print,' '
print,'press enter'
rrr=get_kbrd()
window,3,title='tabulated profile'
plot,relwvl_o,avgobsprof4calib*calcoef1,xtitle='!7Dk!3 [!z(c5)]',ytitle='!6I!3[erg/cm!u2!n/s/sr/Hz]'
oplot,relwvl_o,tabint_to,color=2
xxx_prev=-1
nxwp001z:
print,' '
print,'click left button to mark the wvl point ',iwp+2L
print,'in profile'
print,'press right mouse button to quit marking'
print,' '
lggg_fsac01:
cursor,xxx,yyy,/wait,/data
if(!mouse.button eq 4) then goto,no_nxwp001z
if(xxx eq xxx_prev) then goto,lggg_fsac01 else xxx_prev=xxx
oplot,[xxx,xxx],[0.,1.],color=3,line=2
iwp=iwp+1
rrr=min(abs(xxx-relwvl_o),iwcp)
calibratios[0,iwp]=avgobsprof4calib[iwcp]
calibratios[1,iwp]=tabint_to[iwcp]
calibratios[2,iwp]=xxx+wvlcen_o
calibratios[3,iwp]=iwcp
goto,nxwp001z
no_nxwp001z:
nwp=iwp
calibratios=calibratios[*,0:nwp-1]
calibcoeff_ck=linfit(calibratios[0,*],calibratios[1,*],/double)
window,2
plot,calibratios[0,*],calibratios[1,*],psym=1,$
title='!3calibration curve',xtitle='counts',$
ytitle='erg/cm!u2!n/s/sr/Hz'
oplot,calibratios[0,*],calibratios[0,*]*calibcoeff_ck[1]+calibcoeff_ck[0],color=2
;
;
print,' '
print,' '
print,' '
print,'======================================================'
print,' ABSOLUTE CALIBRATION INTO erg/cm^2/s/sr/Hz '
print,'======================================================'
print,' '
print,'linear calibration coefficients:'
print,'k=',calibcoeff_ck[1]
print,'c=',calibcoeff_ck[0]
scattered_light_disk=fix(((-1.)*calibcoeff_ck[0]/calibcoeff_ck[1])>0.)
print,' '
print,'scattered light at the disc: ',scattered_light_disk,' counts'
print,' '
print,'subtracting the scattered light from the disk area of the image ...'
print,' '
spflfld[*,y1_spdisk:y2_spdisk]=spflfld[*,y1_spdisk:y2_spdisk]-double(scattered_light_disk)
calibratios[0,*]=((calibratios[0,*]-double(scattered_light_disk))>0.)
;
tmpmat1=total(spflfld[twocontareas[0]:twocontareas[1],ypos_sctlofflimb[0]:ypos_sctlofflimb[1]])
n_tmpmat1=double(n_elements(tmpmat1))
tmpmat2=total(spflfld[twocontareas[2]:twocontareas[3],ypos_sctlofflimb[0]:ypos_sctlofflimb[1]])
n_tmpmat2=double(n_elements(tmpmat2))
avg_sctlght_offlimb=(tmpmat1*n_tmpmat1+tmpmat2*n_tmpmat2)/(tmpmat1+tmpmat2)
print,'avg off-limb scattered light: ',fix(avg_sctlght_offlimb),' counts'
print,' '
print,'subtracting the off-limb scattered light ...'
for iy=0L,ydim-1L do begin
if((iy gt y2_spdisk) or (iy lt y1_spdisk)) then begin
spflfld[*,iy]=spflfld[*,iy]-avg_sctlght_offlimb
endif
endfor
spflfld=(spflfld>0.)
loadct,0
window,1,xs=xdim,ys=ydim,title='flatfielded raw spectra -- real relative intensities'
tvscl,spflfld
t_zdg=1.0d0
print,' '
print,'Is a filter put at the disc to decrease its brightness'
print,' to increase contrast of the prominence (y/N)'
sw_sffil=' '
read,sw_sffil
sw_sffil=strupcase(sw_sffil)
if(sw_sffil eq 'Y') then begin
answ_dskftflt=' '
print,' '
print,'Does the filter position fit the disk exactly? (Y/n)'
print,'HINT: If you answer Y the y-positions of the filtered'
print,' area in the spectra image are the same as those'
print,' of the disk area'
read,answ_dskftflt
answ_dskftflt=strupcase(answ_dskftflt)
if(answ_dskftflt ne 'N') then begin
y1_zdgflt=y1_spdisk
y2_zdgflt=y2_spdisk
endif else begin
window,1,xs=xdim,ys=ydim,title='flatfielded raw spectra -- real relative intensities'
tvscl,spflfld
print,'Clicking in the spectra image fix the borders of the filter Y-positions'
print,' '
y1_zfflt=0L
y2_zfflt=0L
print,'click into image to fix the first border'
cursor,xxx,yyy,/wait,/normal
yyy_prev=yyy
y1_zfflt=fix(yyy*float(ydim))
print,y1_zfflt,' pxl'
print,' '
print,'click into image to fix other border'
lgggl_01z01:
cursor,xxx,yyy,/wait,/normal
if(yyy eq yyy_prev) then goto,lgggl_01z01
y2_zfflt=fix(yyy*float(ydim))
print,y2_zfflt
y1_zdgflt=min([y1_zfflt,y2_zfflt])
y2_zdgflt=max([y1_zfflt,y2_zfflt])
endelse
print,' '
print,' '
print,'Filter transmissivity:'
print,'1 ..... 5%'
print,'2 ..... 12.7% (defaults)'
print,'3 ..... other'
zdgflt=' '
read,zdgflt
case zdgflt of
'1' : t_zdg=0.05d0
'2' : t_zdg=0.127d0
'3' : begin
read,t_zdg,prompt='transmissivity in %: '
t_zdg=t_zdg/1.d2
end
else : t_zdg=0.127d0
endcase
spflfld1=spflfld
spflfld1[*,y1_zdgflt:y2_zdgflt]=spflfld[*,y1_zdgflt:y2_zdgflt]/t_zdg
window,1,xs=xdim,ys=ydim,title='flatfielded raw spectra -- real relative intensities'
tvscl,spflfld1
print,' '
print,'View well both filtered and nonfiltered parts? (Y/n)'
read,answ_advv_sw
answ_advv_sw=strupcase(answ_advv_sw)
if(answ_advv_sw ne 'N') then begin
dcspvw=bytarr(3,xdim,ydim)
dcspvw[0,*,0:y1_zdgflt]=bytscl(hist_filter(spflfld1[*,0:y1_zdgflt],percentage=1.))
dcspvw[0,*,y2_zdgflt:*]=bytscl(hist_filter(spflfld1[*,y2_zdgflt:*],percentage=1.))
dcspvw[2,*,y1_zdgflt:y2_zdgflt]=bytscl(hist_filter(spflfld1[*,y1_zdgflt:y2_zdgflt],percentage=1.))
tv,dcspvw,true=1
endif
print,' '
print,'Are you satisfied? (Y/n)'
answ_stsf=' '
read,answ_stsf
answ_stsf=strupcase(answ_stsf)
if(answ_stsf ne 'N') then begin
spflfld=temporary(spflfld1)
endif else begin
spflfld=spflfld_counts
goto,zdgflt_l001
endelse
endif
print,'press enter'
rrr=get_kbrd()
tvlct,rcol2*255,gcol2*255,bcol2*255
print,' '
avgobsprof4calib=total(spflfld[*,ylsp_4calib:yusp_4calib],2)/$
double(yusp_4calib-ylsp_4calib+1L)
for iwp=0,nwp-1L do begin
iwcp=calibratios[3,iwp]
calibratios[0,iwp]=avgobsprof4calib[iwcp]
endfor
sumxy=0.0d0
sumxx=0.0d0
for iwp=0,nwp-1L do begin
sumxy=sumxy+calibratios[0,iwp]*calibratios[1,iwp]
sumxx=sumxx+calibratios[0,iwp]*calibratios[0,iwp]
endfor
calibcoeff_abs=sumxy/sumxx
print,'calibration coefficient after subtraction of the scattered light:'
print,calibcoeff_abs
wset,2
plot,calibratios[0,*],calibratios[1,*],psym=1,$
title='!3calibration curve after scattered-light subtraction',$
xtitle='counts',ytitle='erg/cm!u2!n/s/sr/Hz'
oplot,calibratios[0,*],calibratios[0,*]*calibcoeff_abs,color=2
print,' '
print,' '
print,'Take into account variations of the calib. coeff. with spectral pxls? (Y/n)'
answ_vspxls=' '
read,answ_vspxls
answ_vspxls=strupcase(answ_vspxls)
if(answ_vspxls ne 'N') then begin
;
; variations of calibration coeficient with position on the detector (spectral pixel)
print,' '
print,' '
print,' '
print,'========================================================='
print,' variations of calibration coeficients with position '
print,' on the detector (spectral pixel) '
print,'========================================================='
print,' '
ccfvsppix=reform(calibratios[1,*]/(calibratios[0,*]*calibcoeff_abs))
iscc=where((calibratios[3,*] lt linearea[0]) or (calibratios[3,*] gt linearea[1])) ; <--- exclude line from the fitting
if(iscc[0] ne (-1)) then begin
c_wvlpxlcoeffs=poly_fit(calibratios[3,iscc],ccfvsppix[iscc],2,sigma=errcoeffs,$
chisq=chiwvlpxlfit)
endif else begin
c_wvlpxlcoeffs=poly_fit(calibratios[3,*],ccfvsppix,2,sigma=errcoeffs,$
chisq=chiwvlpxlfit)
endelse
window,8,title='variations of calibration coeficients with position on the detector'
plot,calibratios[2,*],ccfvsppix,xtitle='!7k!3 [!z(c5)]',$
ytitle='!6I!3!dDavid!n/!6I!3!dcalib!n',psym=1
ccffit=c_wvlpxlcoeffs[2]*reform(calibratios[3,*])^2+c_wvlpxlcoeffs[1]*reform(calibratios[3,*])+c_wvlpxlcoeffs[0]
oplot,calibratios[2,*],ccffit,color=2
c_wvlpxlcoeffs=reform(c_wvlpxlcoeffs)
print,' '
print,'a*xpos^2+b*xpos+c'
print,'a=',c_wvlpxlcoeffs[2],' -/+',abs(errcoeffs[2]/c_wvlpxlcoeffs[2]*100.),' %'
print,'b=',c_wvlpxlcoeffs[1],' -/+',abs(errcoeffs[1]/c_wvlpxlcoeffs[1]*100.),' %'
print,'c=',c_wvlpxlcoeffs[0],' -/+',abs(errcoeffs[0]/c_wvlpxlcoeffs[0]*100.),' %'
print,'chi^2=',chiwvlpxlfit
print,' '
print,' '
print,'Are you satisfied? (y/N)'
answ_vspxls1a=' '
read,answ_vspxls1a
answ_vspxls1a=strupcase(answ_vspxls1a)
if(answ_vspxls1a ne 'Y') then begin
if(iscc[0] ne (-1)) then begin
c_wvlpxlcoeffs2=poly_fit(calibratios[3,iscc],ccfvsppix[iscc],1,sigma=errcoeffs,$
chisq=chiwvlpxlfit)
endif else begin
c_wvlpxlcoeffs2=poly_fit(calibratios[3,*],ccfvsppix,1,sigma=errcoeffs,$
chisq=chiwvlpxlfit)
endelse
window,8,title='variations of calibration coeficients with position on the detector'
plot,calibratios[2,*],ccfvsppix,xtitle='!7k!3 [!z(c5)]',$
ytitle='!6I!3!dDavid!n/!6I!3!dcalib!n',psym=1
ccffit=c_wvlpxlcoeffs2[1]*reform(calibratios[3,*])+c_wvlpxlcoeffs2[0]
oplot,calibratios[2,*],ccffit,color=2
c_wvlpxlcoeffs[2]=0.d0
c_wvlpxlcoeffs[1]=c_wvlpxlcoeffs2[1]
c_wvlpxlcoeffs[0]=c_wvlpxlcoeffs2[0]
print,' '
print,'a*xpos+b'
print,'a=',c_wvlpxlcoeffs2[1],' -/+',abs(errcoeffs[1]/c_wvlpxlcoeffs[1]*100.),' %'
print,'b=',c_wvlpxlcoeffs2[0],' -/+',abs(errcoeffs[0]/c_wvlpxlcoeffs[0]*100.),' %'
print,'chi^2=',chiwvlpxlfit
rmvar,c_wvlpxlcoeffs2
print,' '
print,' '
print,'Are you satisfied? (y/N)'
answ_vspxls1a=' '
read,answ_vspxls1a
answ_vspxls1a=strupcase(answ_vspxls1a)
if(answ_vspxls1a ne 'Y') then begin
print,'Then no variations of the calib. coefficient are taken.'
answ_vspxls='N'
endif
endif
endif
if(answ_vspxls eq 'N') then begin
c_wvlpxlcoeffs=dblarr(3)
c_wvlpxlcoeffs[2]=0.
c_wvlpxlcoeffs[1]=0.
c_wvlpxlcoeffs[0]=1.
endif
c_wvlpxlcoeffs=c_wvlpxlcoeffs*calibcoeff_abs
endif
;
; for testing whether the final calibrated profile from postion on the disk
; selected for absolute calibration fit well with David's profile at
; continuum.
;
; tmpvec=spflfld[*,ylsp_4calib:yusp_4calib]
; testprof=total(tmpvec,2)/double(n_elements(tmpvec[0,*]))*(c_wvlpxlcoeffs[0]+c_wvlpxlcoeffs[1]*xpp+$
; c_wvlpxlcoeffs[2]*xpp^2)
; rmvar,tmpvec
; davids=tabint
; iprof=testprof*(c_wvlpxlcoeffs[0]+c_wvlpxlcoeffs[1]*xpp+c_wvlpxlcoeffs[2]*xpp^2)
; save,davids,iprof,wvl,ylsp_4calib,yusp_4calib,file='TEST.idl'
;
;
;
; ************* SAVING THE REDUCED DATA *****************
loadct,0
;
print,' '
print,'Save the flatfielded spectrum into a file? (y/N)'
answ_sffsp=' '
read,answ_sffsp
answ_sffsp=strupcase(answ_sffsp)
if(answ_sffsp eq 'Y') then begin
calib_info='not calibrated'
if(sw_abscal eq 'Y') then begin
calib_info='calibrated into erg/cm^2/s/sr/Hz'
dxpxl_arcsec=one_pxl_in_arcsec ; x&y dimensions of one pixel in arcsec
dypxl_arcsec=one_pxl_in_arcsec
xcent_pxl=disccenter_pxl[0] ; X&Y coordinates of the disk center in pxl
ycent_pxl=disccenter_pxl[1]
uwire_sj_pxl=ypos_wires_sj[*,1] ; [x,y] positions in pxl of intersections of the upper horizontal wire with the slit
lwire_sj_pxl=ypos_wires_sj[*,0] ; [x,y] positions in pxl of intersections of the upper horizontal wire with the slit
endif
print,' '
fltype=' '
path_out=path_sp
print,' '
print,'The file(s) are gonna be saved into the directory:'
print,path_out
print,' '
print,'Do you wanna change this directory? (y/N)'
answ_choutdir=' '
read,answ_choutdir
answ_choutdir=strupcase(answ_choutdir)
if(answ_choutdir eq 'Y') then path_out=dialog_pickfile(path='c:\MFS_zpracovani\',/directory,title='directory where to save file(s)')
plnth=strlen(path_out)
if(strmid(path_out,plnth-1L) ne dir_separator) then path_out=path_out+dir_separator
if(sw_wvl eq 'Y') then begin
answ_z001=' '
print,'Save as:'
print,'1 .... idl-save file(s) (default)'
print,'2 .... fits file(s) with standard header'
read,answ_z001
ex_answ_z001=where(answ_z001 eq ['1','2'])
if(ex_answ_z001[0] eq (-1)) then answ_z001='1'
if(answ_z001 eq '1') then fltype='idl'
if(answ_z001 eq '2') then fltype='fts'
endif else begin
fltype='idl'
endelse
;
date_sp=sd3+monnum2str(fix(sd2))+sd1a
time_sp=st1+':'+st2+':'+st3+'.'+st4
if(fltype eq 'idl') then begin
if(sw_wvl eq 'Y') then begin
if(sw_abscal eq 'Y') then begin
flspout=path_out+flnm_sp+'_ff_abscalib.idl'
spflfld_calibrated=dblarr(xdim,ydim)
xcoords_sj_arcsec=(findgen(xdim_sj)-xcent_pxl)*dxpxl_arcsec
ycoords_sj_arcsec=(findgen(ydim_sj)-ycent_pxl)*dypxl_arcsec
for iy=0L,ydim-1L do begin
tmpvec1=reform(spflfld[*,iy])*(c_wvlpxlcoeffs[0]+c_wvlpxlcoeffs[1]*xpp+c_wvlpxlcoeffs[2]*xpp^2)
spflfld_calibrated[*,iy]=tmpvec1
endfor
save,spflfld_counts,spflfld_calibrated,wvl,date_sp,time_sp,ypos_wires_sp,calib_info,imgsj,$
xcoords_sj_arcsec,ycoords_sj_arcsec,date_sj,time_sj,uwire_sj_pxl,$
lwire_sj_pxl,file=flspout
endif else begin
flspout=path_out+flnm_sp+'_ff_wvlcalib.idl'
save,spflfld,wvl,date_sp,time_sp,ypos_wires_sp,calib_info,file=flspout
endelse
endif else begin
flspout=path_out+flnm_sp+'_ff.idl'
save,spflfld,date_sp,time_sp,ypos_wires_sp,calib_info,file=flspout
endelse
print,' '
print,'% the file '
print,'% '+flspout+' saved'
print,' '
endif
if(fltype eq 'fts') then begin
if(sw_abscal eq 'Y') then flspout=path_out+flnm_sp+'_HalpSPint.fts' else $
flspout=path_out+flnm_sp+'_ff_wvlcalib.fts'
outhdr=sphdr
fits_open,flspout,fcbout,/write
sxaddpar,outhdr,'DATE-OBS',date_sp,''
sxaddpar,outhdr,'TIME',time_sp,''
sxaddpar,outhdr,'UPHAIR',ypos_wires_sp[1],'y position of top hair on the image [pxl]'
sxaddpar,outhdr,'DOWNHAIR',ypos_wires_sp[0],'y position of bottom hair on the image [pxl]'
sxaddpar,outhdr,'LAMBDA 0',round(pxl2wvl[0]*1.e3)/1.e3,'wavelength of the first pixel [Angstroems]'
sxaddpar,outhdr,'DISPCOEF',round(pxl2wvl[1]*1.e4)/1.e4,'dispersion curve -- linear coeff [Angstroems/pxl]'
if(sw_abscal eq 'Y') then begin
sxaddpar,outhdr,'INT0COEF',c_wvlpxlcoeffs[0],'calib coef into erg/cm^2/s/sr/Hz of the 1st pxl'
sxaddpar,outhdr,'INT1COEF',c_wvlpxlcoeffs[1],'calib coef - linear change with spectral pixels'
sxaddpar,outhdr,'INT2COEF',c_wvlpxlcoeffs[2],'calib coef - quadratic change with spectral pixels'
sxaddpar,outhdr,'COMMENT','X -- pos. in wvl pxls starting from zero from left edge of the image'
sxaddpar,outhdr,'COMMENT','Y -- pos. in pxls along the slit starting from zero from the bottom'
sxaddpar,outhdr,'COMMENT','val[X,Y] -- value at pixel position [X,Y]'
sxaddpar,outhdr,'COMMENT','wavelength in Angstroems: X*DISPCOEF+LAMBDA 0'
sxaddpar,outhdr,'COMMENT','calibrated intensity in erg/cm^2/s/sr/Hz'
sxaddpar,outhdr,'COMMENT','val[X,Y]*(INT0COEF+INT1COEF*X+INT2COEF*x^2)'
endif
sxaddpar,outhdr,'CALIBINF',calib_info,' '
fits_write,fcbout,spflfld,outhdr
fits_close,fcbout
print,' '
print,'% the file '
print,'% '+flspout+' saved'
print,' '
if(sw_abscal eq 'Y') then begin
flspout_counts=path_out+flnm_sp+'_HalpSPcnt.fts'
outhdr=sphdr
fits_open,flspout_counts,fcbout,/write
sxaddpar,outhdr,'DATE-OBS',date_sp,''
sxaddpar,outhdr,'TIME',time_sp,''
sxaddpar,outhdr,'UPHAIR',ypos_wires_sp[1],'y position of top hair on the image [pxl]'
sxaddpar,outhdr,'DOWNHAIR',ypos_wires_sp[0],'y position of bottom hair on the image [pxl]'
sxaddpar,outhdr,'LAMBDA 0',round(pxl2wvl[0]*1.e3)/1.e3,'wavelength of the first pixel [Angstroems]'
sxaddpar,outhdr,'DISPCOEF',round(pxl2wvl[1]*1.e4)/1.e4,'dispersion curve -- linear coeff [Angstroems/pxl]'
sxaddpar,outhdr,'COMMENT','X -- pos. in wvl pxls starting from zero from left edge of the image'
sxaddpar,outhdr,'COMMENT','Y -- pos. in pxls along the slit starting from zero from the bottom'
sxaddpar,outhdr,'COMMENT','val[X,Y] -- observed intensity in counts at pixel position [X,Y]'
sxaddpar,outhdr,'COMMENT','wavelength in Angstroems: X*DISPCOEF+LAMBDA 0'
fits_write,fcbout,spflfld_counts,outhdr
fits_close,fcbout
print,' '
print,'% the file '
print,'% '+flspout_counts+' saved'
print,' '
;
flsjout=path_out+flnm_sp+'_HalpSJ.fts'
outhdr1=hdrsj
fits_open,flsjout,fcbout,/write
sxdelpar,outhdr1,'DATE'
sxaddpar,outhdr1,'DATE-OBS',date_sj,''
sxaddpar,outhdr1,'TIME',time_sj,''
sxaddpar,outhdr1,'UPHAIRX',fix(round(uwire_sj_pxl[0])),'x position of top hair on the slit [pxl]'
sxaddpar,outhdr1,'UPHAIRY',fix(round(uwire_sj_pxl[1])),'y position of top hair on the slit [pxl]'
sxaddpar,outhdr1,'DWNHAIRX',fix(round(lwire_sj_pxl[0])),'x position of bottom hair on the slit [pxl]'
sxaddpar,outhdr1,'DWNHAIRY',fix(round(lwire_sj_pxl[1])),'y position of bottom hair on the slit [pxl]'
sxaddpar,outhdr1,'CRPIX1',fix(round(xcent_pxl)),'x position of the disk center [pxl]'
sxaddpar,outhdr1,'CRPIX2',fix(round(ycent_pxl)),'y position of the disk center [pxl]'
sxaddpar,outhdr1,'CDELT1',round(dxpxl_arcsec*1.e4)/1.e4,'x-size of the pixel [arcsec]'
sxaddpar,outhdr1,'CDELT2',round(dypxl_arcsec*1.e4)/1.e4,'y-size of the pixel [arcsec]'
fits_write,fcbout,imgsj,outhdr1
fits_close,fcbout
print,' '
print,'% the file '
print,'% '+flsjout+' saved'
print,' '
endif
endif
endif
endif
;
eocmffffs:
end