; ; ; ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ; + + ; + 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)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_s10.)<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