;
;
;   ++++++++++++++++++++++++++++++++++++++++++++++++++++++++
;   +                                                      +
;   +               FUNCTIONS & PROCEDURES                 +
;   +                                                      +
;   ++++++++++++++++++++++++++++++++++++++++++++++++++++++++
;
;
pro make_setup
;
common setup_cb,idl_4degpolyfitting_sw,polyfitting_path,basepath
;
; 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
;
; basepath .... path to all file needed by the code; at the end dir sperator is obligatory, thus path should 
;                look like this, e.g., 'c:\MFS_zpracovani\MFS_zpracovani_mg\new\'
basepath='c:\MFS_zpracovani\uni_preparation\'
;
; 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=basepath ; <--- 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)<top_limit)>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)<ul)
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<mm[1])
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<mm[1])
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 calibrefprof,tabfile,mu,wvl,int,relwvl=relwvl,symmetric_from_halfprof=symmetric_from_halfprof,problems=problems
;
problems=''
mu=double(mu)
print,'***: ',mu
if(mu gt 1.) then begin
 print,'%'
 print,'% calibrefprof: mu cannot be gt 1'
 print,'%               value of mu set to 1.0'
 print,'%'
 mu1=1.d0
 problems=problems+' mu too high'
endif else begin
 if(mu lt 0.141) then begin
  print,'%'
  print,'% calibrefprof: to close to the limb. Value of mu'
  print,'%               cahged to 0.141. '
  mu1=.141d0
  problems=problems+' mu too low'
 endif else begin
  mu1=mu
 endelse
endelse
centwvl=0.
imax=0.
nnn=0
;
if(n_elements(symmetric_from_halfprof) eq 0) then symmetric_from_halfprof=0
;values of mu for 7 columns of intensities
musforcols=[1.0d0,0.8d0,0.6d0,0.436d0,0.312d0,0.28d0,0.141d0]
;
openr,1,tabfile
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,*)
for i=0,nnn-1 do begin
 lin_interp,musforcols,reform(intens(*,i)),[mu1],profval
 halfprof(1,i)=profval
endfor
if(symmetric_from_halfprof gt 0) then begin
 zzz=where(halfprof eq 0)
 if(zzz[0] ne (-1)) then begin
  complprof=fltarr(2,nnn*2L-1L)
  complprof(0,0:nnn-1L)=reverse(halfprof(0,*),2)*(-1.0)
  complprof(1,0:nnn-1L)=reverse(halfprof(1,*),2) 
  complprof(*,nnn:nnn*2L-2L)=halfprof(*,1:nnn-1L)
  complprof(0,*)=complprof(0,*)+centwvl
 endif else begin  
  complprof=fltarr(2,nnn*2L)
  complprof(0,0:nnn-1L)=reverse(halfprof(0,*),2)*(-1.0)
  complprof(1,0:nnn-1L)=reverse(halfprof(1,*),2) 
  complprof(*,nnn:nnn*2L-1L)=halfprof(*,0:nnn-1L)
  complprof(0,*)=complprof(0,*)+centwvl
 endelse
endif else begin
 complprof=fltarr(2,nnn)
 complprof(0,*)=halfprof(0,*)+centwvl
 complprof(1,*)=halfprof(1,*)
endelse
lin_interp,reform(complprof(0,*)),reform(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=-1L
i_abs=-1L
;
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
   i_emis=i_emis+1L
   minemisvec[i_emis]=avgcval
  endif
  if(avglval lt avgcval) then begin ; absorption profile
   i_abs=i_abs+1
   maxabsvec[i_abs]=avgcval
  endif
 endif
endfor
minemis=0.0d0
maxabs=0.0d0
ggg=0
if(i_emis ge 0) then begin
 minemisvec=minemisvec[0:i_emis]
 minemis=min(minemisvec)
endif else begin
 print,'absorption profiles only !!!!'
 ggg=-1
endelse
if(i_abs ge 0) then begin
 maxabsvec=maxabsvec[0:i_abs]
 maxabs=max(maxabsvec)
endif else begin
 ggg=1
 print,'emission profiles only !!!!'
endelse
;
;
zerswsw=0
if(minemis le 0.) then zerswsw=zerswsw+1
if(maxabs le 0.) then zerswsw=zerswsw+1
if(zerswsw eq 0) then begin
 if(ggg eq 0) then ratio_maxabs_minemis=minemis/maxabs else ratio_maxabs_minemis=0.0
endif else begin
 print,'negative values in some profiles!!!!!'
 ratio_maxabs_minemis=0.0
endelse
;
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
;
;
;
pro load_calibparams,file_in,c0,c1,c2,wvl0,disp,sctlght_disk=sctlght_disk,$
 sctlght_offlimb=sctlght_offlimb,twocontsections=twocontsections
;
fits_read,file_in,data,header
c0=fxpar(header,'INT0COEF')
c1=fxpar(header,'INT1COEF')
c2=fxpar(header,'INT2COEF')
wvl0=fxpar(header,'LAMBDA 0')
disp=fxpar(header,'DISPCOEF')
sctlght_disk=fxpar(header,'SCTLDSK')
sctlght_offlimb=fxpar(header,'SCTLLMB')
twocontsections=lonarr(4)
twocontsections[0]=fxpar(header,'CNTARLL')
twocontsections[1]=fxpar(header,'CNTARLR')
twocontsections[2]=fxpar(header,'CNTARRL')
twocontsections[3]=fxpar(header,'CNTARRR')
;
end
;
;
;
pro filename2dandt,flnm,datestr,timestr
;
posd1=strpos(flnm,'T',/reverse_search)
if(posd1[0] eq (-1)) then posd1=strpos(flnm,'_',/reverse_search)
sd1=strmid(flnm,posd1-6,2)
sd2=strmid(flnm,posd1-4,2)
sd3=strmid(flnm,posd1-2,2)
st1=strmid(flnm,posd1+1,2)
st2=strmid(flnm,posd1+3,2)
st3=strmid(flnm,posd1+5,2)
tst1=strmid(flnm,posd1+7,1)
tst1=strupcase(tst1)
tst1t=where(tst1 eq ['N','F','D','.','_'])
if(tst1t[0] eq (-1)) then begin
 st4=strmid(flnm,posd1+7,2)
endif else begin
 st4='00'
endelse
if(fix(sd1) gt 50) then sd1a='19'+sd1 else sd1a='20'+sd1
datestr=' '
timestr=' '
datestr=sd3+'-'+monnum2str(fix(sd2))+'-'+sd1a
timestr=st1+':'+st2+':'+st3+'.'+st4
end
;
;
;
;
;
;
;         +++++++++++++++++++++++++++++++++++++++++
;         +                                       +
;         +               MAIN CODE               +
;         +                                       +
;         +++++++++++++++++++++++++++++++++++++++++
;
;
common setup_cb,idl_4degpolyfitting_sw,polyfitting_path,basepath
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
device,get_screen_size=wsspdd
maxywinsize=round(float(wsspdd[1])*0.8)
;
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=basepath+'raw_data',/directory)
  fls_dc=file_search(fls_dcd,'*D.fts',count=n_fls_dc)
 endif else begin
  fls_dc=dialog_pickfile(path=basepath+'raw_data',/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=basepath+'raw_data',/directory)
  fls_ff=file_search(fls_ffd,'*F.fts',count=n_fls_ff)
 endif else begin
  fls_ff=dialog_pickfile(path=basepath+'raw_data',/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)
 wxdim=800
 wydim=round(float(800)*float(ydim)/float(xdim))
 if(wydim gt maxywinsize) then begin
  rrrw=float(xdim)/float(ydim)
  wydim=maxywinsize
  wxdim=round(float(wydim)*rrrw)
 endif
 window,0,xs=wxdim,ys=wydim,title='mean FLAT FRAME'
 showimg,data=congrid(ffmat,wxdim,wydim)
 xloadct,/block, UPDATECALLBACK='showimg',UPDATECBDATA=congrid(ffmat,wxdim,wydim)
 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=''
 read,rrr,prompt=''
 window,0,xs=wxdim,ys=wydim,title='mean FLAT FRAME'
 showimg,data=congrid(ffmat,wxdim,wydim)
 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=' '
 read,rrr,prompt=''
 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=wxdim,ys=wydim,title='mean FLAT FRAME'
 tv,congrid(ffmatb,wxdim,wydim)
 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=' '
 read,rrr,prompt=''
 window,0,xs=wxdim,ys=wydim,title='mean FLAT FRAME'
 tv,congrid(ffmatb,wxdim,wydim)
 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=wxdim,ys=wydim,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,congrid(ffmatb,wxdim,wydim)
 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<mm[1])
 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=wxdim,ys=wydim,title='mean FLAT FRAME corrected for inclination'
 tvscl,congrid(ffmat_s2,wxdim,wydim)
 ;
 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=' '
 read,rrr,prompt=''
 window,0,xs=wxdim,ys=wydim,title='mean FLAT FRAME corrected for inclination'
 tvscl,congrid(ffmat_s2,wxdim,wydim)
 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=wxdim,ys=wydim,title='mean FLAT FRAME corrected for both inclination & curvature'
 tvscl,congrid(ffmat_s3,wxdim,wydim)
 wshow,0
 rrr=' '
 read,rrr,prompt=''
 ;
 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=wxdim,ys=wydim,title='soft flat matrix'
 tvscl,congrid(hist_filter(ffmat_ms,percentage=1.),wxdim,wydim)
 ;
 ;
 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=wxdim,ys=wydim,title='hard flat matrix'
 tvscl,congrid(dtfilter(ffmat_mh,threshold=2.),wxdim,wydim)
 print,' '
 print,'press enter'
 rrr=' '
 read,rrr,prompt=''
 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=wxdim,ys=wydim,title='slit-flat matrix (drifting)'
 tvscl,congrid(dtfilter(ffmat_slit,threshold=2.),wxdim,wydim)
 wset,1
 tvscl,congrid(dtfilter(ffmat_mhp,threshold=2.),wxdim,wydim)
 print,' '
 print,'press enter'
 rrr=' '
 read,rrr,prompt=''
 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=wxdim,ys=wydim,title='static-flat matrix'
 tvscl,congrid(hist_filter(ffmat_static,percentage=1.),wxdim,wydim)
 ;
 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=basepath+'raw_data',/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=basepath+'calibrated',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=basepath+'raw_data',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)
 wxdim=800
 wydim=round(float(800)*float(ydim)/float(xdim))
 if(wydim gt maxywinsize) then begin
  rrrw=float(xdim)/float(ydim)
  wydim=maxywinsize
  wxdim=round(float(wydim)*rrrw)
 endif
 window,0,xs=wxdim,ys=wydim,title='raw spectra'
 showimg,data=congrid(spmat,wxdim,wydim)
 xloadct,/block, UPDATECALLBACK='showimg',UPDATECBDATA=congrid(spmat,wxdim,wydim)
 tvlct,rcol_sp,gcol_sp,bcol_sp,/get
 ; horizontal wires
 print,' '
 print,'Click on the two horizontal spectrograph wires'
 print,'press enter'
 rrr=' '
 read,rrr,prompt=''
 window,0,xs=wxdim,ys=wydim,title='raw spectra'
 showimg,data=congrid(spmat,wxdim,wydim)
 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=' '
 read,rrr,prompt=''
 ypos_wires_sp=[min([ypw1,ypw2]),max([ypw1,ypw2])]
 print,'% subtracting dark-frame'
 spmat1=spmat-dcmat
 n_below_zero=where(spmat1 lt 0.d0)
 if(n_below_zero[0] ne (-1)) then begin
  n_perc_below_zero=double(n_elements(n_below_zero))/double(xdim*ydim)*100.d0
  spmat1=(spmat1>0.)
  print,' '
  if(n_perc_below_zero lt 30.d0) then begin
   if(n_perc_below_zero lt 0.1d0) then begin
    print,'After subtraction of avg dark-frame there is less then 0.1% of pxls with'
    print,' count values below zero. A zero value was set for these pxls.'
    wait,0.2
   endif else begin
    print,'After subtraction of avg dark-frame there is '+string(n_perc_below_zero,format='(f5.1)')+'% of pxls with'
    print,' count values below zero. A zero value was set for these pxls.'
    wait,0.2
   endelse
  endif else begin
  ;tuto
   print,'WARNING: After subtraction of avg dark-frame there is '+string(n_perc_below_zero,format='(f5.1)')+'% of pxls with'
   print,'         count values below zero. It means that the dark current is too high for'
   print,'         these observations. A zero value was set for these pxls.'
   print,' '
   answ_exxf=''
   print,'Exit or Continue? (X or C)' 
   l_eorc_ac2:
   read,answ_exxf
   answ_exxf=strcompress(answ_exxf,/remove_all)
   answ_exxf=strmid(answ_exxf,0,1)
   answ_exxf=strupcase(answ_exxf)
   rrr=where(answ_exxf eq ['X','C'])
   if(rrr[0] eq (-1)) then begin
    print,' '
    print,'Please, answer ''c'' or ''x''!!!!'
    print,' '
    goto,l_eorc_ac2
   endif
   if(answ_exxf eq 'X') then begin 
    print,'Exiting ...'
    goto,eocmffffs
   endif
  endelse
  print,' '
 endif
 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.)
 window,0,xs=wxdim,ys=wydim,title='corrected raw spectra'
 showimg,data=congrid(spmat3,wxdim,wydim)
 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
 if(ffydrift ne 0.d0) then begin
  print,' '
  print,'Apply the drift? (y/N)'
  ydrft_sw=' '
  read, ydrft_sw
  ydrft_sw=strcompress(ydrft_sw,/remove_all)
  ydrft_sw=strmid(ydrft_sw,0,1)
  ydrft_sw=strupcase(ydrft_sw)
 endif else begin
  ydrft_sw='N'
 endelse
 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.
 fltfield=dtfilter(fltfield,threshold=4.)
 spflfld=spmat3/fltfield ; flatfielded spectrum
 tvlct,rcol_sp,gcol_sp,bcol_sp
 window,0,xs=wxdim,ys=wydim,title='corrected raw spectra'
 tvscl,congrid(spmat3,wxdim,wydim) 
 window,1,xs=wxdim,ys=wydim,title='flatfielded raw spectra'
 tvscl,congrid(spflfld,wxdim,wydim)
 ad_vis_adj_sw=''
 print,' '
 print,'Ajust visibility of the flatfielded spectra? (y/N)'
 read,ad_vis_adj_sw
 ad_vis_adj_sw=strcompress(ad_vis_adj_sw,/remove_all)
 ad_vis_adj_sw=strmid(ad_vis_adj_sw,0,1)
 ad_vis_adj_sw=strupcase(ad_vis_adj_sw)
 if(ad_vis_adj_sw eq 'Y') then begin 
  xloadct,/block, UPDATECALLBACK='showimg',UPDATECBDATA=congrid(spflfld,wxdim,wydim)
  tvlct,rcol_sp,gcol_sp,bcol_sp,/get
 endif
 loadct,0
 window,2,xs=wxdim,ys=wydim,title='flatfield'
 tvscl,congrid(hist_filter(fltfield),wxdim,wydim)
 ;
 print,' '
 print,' '
 print,'For calibration it is necessary that the spectrograph slit'
 print,' intersects the disk. If this is not the case, parameters'
 print,' of wvl- and intensity-calibration curves should be taken'
 print,' from fits-file with calibrated disk-data.'
 print,' '
 print,'Does the slit intersect the solar disk? (Y/n)'
 sltintsc_dsk=' '
 read,sltintsc_dsk
 print,' '
 print,' '
 sltintsc_dsk=strupcase(sltintsc_dsk)
 tmpval=where(sltintsc_dsk eq ['Y','N'])
 if(tmpval[0] eq (-1)) then sltintsc_dsk='Y'
 wdelete,0
 wdelete,2
 if(sltintsc_dsk eq 'N') then begin 
  fl_dskcalib=dialog_pickfile(title='Fits-file with calibrated disk-data')
  load_calibparams,fl_dskcalib,c0,c1,c2,wvl_o,dispersion,$
                 sctlght_disk=scattered_light_disk,sctlght_offlimb=avg_sctlght_offlimb,$
                 twocontsections=twocontareas
  print,'If you had to reverse in wvl the spectrum in disk-data'
  print,' the prominence spectrum should be reversed as well.'
  print,'Reverse the spectrum in wvl? (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)
   tvlct,rcol_sp,gcol_sp,bcol_sp
   window,1,xs=wxdim,ys=wydim,title='flatfielded raw spectra'
   tvscl,congrid(spflfld1,wxdim,wydim)
   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
   print,' '
  endif
  sep_pathflnm,fl_dskcalib,path_dskcalib,flnm_dskcalib,flext_dskcalib
  print,' '
  print,'++++++++++++++++++++++++++++++++++++++++++++++++++++++++'
  print,'coefficients for wvl calibration retreaved from file '
  print,flnm_dskcalib+'.'+flext_dskcalib+' are: '
  print,'wvl_o: ',wvl_o
  print,'dispersion: ',dispersion
  pxl2wvl=dblarr(2)
  pxl2wvl[0]=wvl_o
  pxl2wvl[1]=dispersion
  print,' '
  print,'intensity calibration'
  print,'C_o: ',c0
  print,'C_1: ',c1
  print,'C_2: ',c2
  c_wvlpxlcoeffs=dblarr(3)
  c_wvlpxlcoeffs[0]=c0
  c_wvlpxlcoeffs[1]=c1
  c_wvlpxlcoeffs[2]=c2
  print,'scattered light at disk: '+numb2string(scattered_light_disk)+' counts'
  print,'off-limb scattered light: '+numb2string(avg_sctlght_offlimb)+' counts'
  print,'++++++++++++++++++++++++++++++++++++++++++++++++++++++++'
  print,' '
  spflfld_counts=spflfld
  answ_estolsctlght=' '
  print,'Estimate avg intensity of off-limb scattered light? (Y/n)'
  print,'If you answere ''n'' value from disk-data file is taken.'
  read,answ_estolsctlght
  answ_estolsctlght=strupcase(answ_estolsctlght)
  if(answ_estolsctlght ne 'N') then begin
   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=' '
   read,rrr,prompt=''
   tvlct,rcol_sp,gcol_sp,bcol_sp
   window,1,xs=wxdim,ys=wydim,title='flatfielded raw spectra'
   tvscl,congrid(spflfld,wxdim,wydim)
   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_z01ab35:
   cursor,xxx,yyy,/wait,/normal
   if(yyy eq yyy_prev) then goto,lggg_z01ab35
   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])
   tmpmat1=total(spflfld[twocontareas[0]:twocontareas[1],*])
   n_tmpmat1=double(n_elements(tmpmat1))
   tmpmat2=total(spflfld[twocontareas[2]:twocontareas[3],*])
   n_tmpmat2=double(n_elements(tmpmat2))
   avg_sctlght_offlimb=(tmpmat1*n_tmpmat1+tmpmat2*n_tmpmat2)/(tmpmat1+tmpmat2)
   avg_sctlght_offlimb=round(avg_sctlght_offlimb)
   print,'avg off-limb scattered light: ',fix(avg_sctlght_offlimb),' counts'
  endif
  print,' '
  print,'subtracting the off-limb scattered light ...'
  spflfld=spflfld-avg_sctlght_offlimb
  tvlct,rcol_sp,gcol_sp,bcol_sp
  window,1,xs=wxdim,ys=wydim,title='flatfielded spectra, scattered light removed'
  tvscl,congrid(spflfld,wxdim,wydim)
  ;
  print,' '
  print,' '
  print,'Pick the slit-jaw complementary to the spectra'
  print,' read from the file '+flnm_sp+'.'+flext_sp
  fl_sj=dialog_pickfile(title='Fits-file with corresp.slit-jaw')
  fits_read,fl_sj,imgsj,hdrsj
  sep_pathflnm,fl_sj,path_sj,flnm_sj,flext_sj
  sss=size(imgsj)
  xdim_sj=sss[1]
  ydim_sj=sss[2]
  wxdim_sj=800
  wydim_sj=round(float(wxdim_sj)*float(ydim_sj)/float(xdim_sj))
  if(wydim_sj gt maxywinsize) then begin
   wydim_sj=maxywinsize
   wxdim_sj=round(float(wydim_sj)*float(xdim_sj)/float(ydim_sj))
  endif
  loadct,0   
  window,4,xs=wxdim_sj,ys=wydim_sj,title='slit-jaw image'
  showimg,data=congrid(imgsj,wxdim_sj,wydim_sj)
  xloadct,/block, UPDATECALLBACK='showimg',UPDATECBDATA=congrid(imgsj,wxdim_sj,wydim_sj)
  tvlct,rcol_sj,gcol_sj,bcol_sj,/get
  ; ---> Maybe, I will insert some other portion of code here
  filename2dandt,flnm_sj,date_sj,time_sj ; tuto
  tvlct,rcol_sj,gcol_sj,bcol_sj
  tvscl,congrid(imgsj,wxdim_sj,wydim_sj)
  contour,imgsj,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'
  window,4,xs=wxdim_sj,ys=wydim_sj,title='slit-jaw image'
  tvscl,congrid(imgsj,wxdim_sj,wydim_sj)
  contour,imgsj,findgen(xdim_sj),findgen(ydim_sj),ticklen=-1,$
   position=[0,0,xdim_sj-1,ydim_sj-1],xstyle=1,ystyle=1,/noerase,/nodata
  lggg_thwsp08az:
  cursor,xx,yy,/wait,/data
  if((xx eq xx_prev) and (yy eq yy_prev)) then goto,lggg_thwsp08az
  print,xx,' pxl   ',yy,' pxl'
  print,' '
  print,' '
  ypos_wires_sj[0,0]=xx
  ypos_wires_sj[1,0]=yy
  if(ypos_wires_sj[1,0] gt ypos_wires_sj[1,1]) then begin
   tmpval01ulws=ypos_wires_sj[1,0]
   ypos_wires_sj[1,0]=ypos_wires_sj[1,1]
   ypos_wires_sj[1,1]=tmpval01ulws
   rmvar,tmpval01ulws
  endif
  xcent_pxl=0
  ycent_pxl=0
  dxpxl_arcsec=0.4507
  dypxl_arcsec=0.4507
  sw_wvl='Y'
  sw_abscal='Y'
  only_disk=-1
  goto,l_svresults
 endif
 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,basepath+'Mg1_catprof.idl'

 window,0,title='sample avg disc spectrum [ window No.1 ]'
 plot,mg1_catprof[0,*],mg1_catprof[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)
  tvlct,rcol_sp,gcol_sp,bcol_sp
  window,1,xs=wxdim,ys=wydim,title='flatfielded raw spectra'
  tvscl,congrid(spflfld1,wxdim,wydim)
  loadct,0
  window,0,title='sample avg disc spectrum'
  plot,mg1_catprof[0,*],mg1_catprof[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'
   tvlct,rcol_sp,gcol_sp,bcol_sp
   window,1,xs=wxdim,ys=wydim,title='flatfielded raw spectra'
   tvscl,congrid(spflfld,wxdim,wydim)
   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,mg1_catprof[0,*],mg1_catprof[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=' '
  read,rrr,prompt=''
  wshow,2
  print,' '
  print,'Now, mark the same wavelength points in the sample spectrum'
  window,0,title='sample avg disc spectrum'
  plot,mg1_catprof[0,*],mg1_catprof[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=' '
  read,rrr,prompt=''
  wdelete,0
 endif
 loadct,0
 print,' '
 print,' '
 ;
 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:
  tvlct,rcol_sp,gcol_sp,bcol_sp
  window,1,xs=wxdim,ys=wydim,title='flatfielded raw spectra'
  tvscl,congrid(spflfld,wxdim,wydim)
  answ_advv_sw='N'
  tmpstr=''
  print,' '
  print,'---------------------------------------------------------'
  print,'1 ........ observation at the limb (both disk and limb)'
  print,'           (default)'
  print,' '
  print,'2 .......  observation at the disk (only disk)'
  print,' '
  read,tmpstr
  print,' '
  tmpstr=strcompress(tmpstr,/remove_all)
  tmpstr=tmpstr[0]
  tmpval=where(tmpstr eq ['1','2'])
  if(tmpval[0] eq (-1)) then tmpstr='1'
  if(tmpstr eq '1') then only_disk=0
  if(tmpstr eq '2') then only_disk=1
  if(only_disk eq 1) then begin
   muest=0.
   print,'mu for observed position on the disk (1 for the disk center)'
   read,muest
   y1_spdisk=0L
   y2_spdisk=ydim-1L
   spflfld_counts=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'
  endif else begin
   print,'please, chose Y borders of the disk in the spectrum image'
   print,' '
   print,'press enter'
   rrr=' '
   read,rrr,prompt=''
   tvlct,rcol_sp,gcol_sp,bcol_sp
   window,1,xs=wxdim,ys=wydim,title='flatfielded raw spectra'
   tvscl,congrid(spflfld,wxdim,wydim)
   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=basepath,title='Pick the slit-jaw fits file')
   fits_read,fl_sj,imgsj,hdrsj
   sep_pathflnm,fl_sj,path_sj,flnm_sj,flext_sj
   sss=size(imgsj)
   xdim_sj=sss[1]
   ydim_sj=sss[2]
   loadct,0
   wxdim_sj=800
   wydim_sj=round(float(wxdim_sj)*float(ydim_sj)/float(xdim_sj))
   if(wydim_sj gt maxywinsize) then begin
    wydim_sj=maxywinsize
    wxdim_sj=round(float(wydim_sj)*float(xdim_sj)/float(ydim_sj))
   endif
   window,4,xs=wxdim_sj,ys=wydim_sj,title='slit-jaw image'
   showimg,data=congrid(imgsj,wxdim_sj,wydim_sj)
   xloadct,/block, UPDATECALLBACK='showimg',UPDATECBDATA=congrid(imgsj,wxdim_sj,wydim_sj)
   tvlct,rcol_sj,gcol_sj,bcol_sj,/get
   rcol_sj[255]=255
   gcol_sj[255]=255
   bcol_sj[255]=255
   rcol_sj[254]=255
   gcol_sj[254]=0
   bcol_sj[254]=0
   rcol_sj[253]=200
   gcol_sj[253]=200
   bcol_sj[253]=0
   tvlct,rcol_sj,gcol_sj,bcol_sj
   imgsjb=bytscl(imgsj,top=252)
   estcc_l001:
   window,4,xs=wxdim_sj,ys=wydim_sj,title='slit-jaw image'
   tv,congrid(imgsjb,wxdim_sj,wydim_sj)
   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
   if(ypos_wires_sj[1,0] gt ypos_wires_sj[1,1]) then begin
    tmpval01ulws=ypos_wires_sj[1,0]
    ypos_wires_sj[1,0]=ypos_wires_sj[1,1]
    ypos_wires_sj[1,1]=tmpval01ulws
    rmvar,tmpval01ulws
   endif
   window,4,xs=wxdim_sj,ys=wydim_sj,title='slit-jaw image'
   tv,congrid(imgsjb,wxdim_sj,wydim_sj)
   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=' '
   read,rrr,prompt=''
   window,4,xs=wxdim_sj,ys=wydim_sj,title='slit-jaw image'
   tv,congrid(imgsjb,wxdim_sj,wydim_sj)
   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,' '
   tvlct,rcol_sp,gcol_sp,bcol_sp
   window,1,xs=wxdim,ys=wydim,title='flatfielded raw spectra -- real relative intensities'
   tvscl,congrid(spflfld,wxdim,wydim)
   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,' '
   filename2dandt,flnm_sj,date_sj,time_sj
   cmpltdtstr=date_sj+' '+time_sj
   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
  endelse
  ;
  if(only_disk eq 0) then begin
   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=''
   read,rrr,prompt=''
   window,1,xs=wxdim,ys=wydim,title='flatfielded raw spectra -- real relative intensities'
   tvscl,congrid(spflfld,wxdim,wydim)
   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])
  endif
  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=' '
  read,rrr,prompt=''
  loadct,0
  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
  calibrefprof,basepath+'MG1_5172_7A_accord2cont.tab',muest,wvl,tabint,relwvl=relwvl_t,problems=clbrp_problems
  rrrcrbbp=strpos(clbrp_problems,'mu too low')
  if(rrrcrbbp ne (-1)) then begin
   print,''
   print,'WARNING:'
   print,'Slit too close to the limb. It is recomended to make calibration of '
   print,' spectra from the disc center. And load calibration coefficients from '
   print,' FITS file of the calibrated spectrum from  the disc center.'
   answ_exxf=''
   print,'Exit or Continue? (X or C)' 
   l_eorc_ac:
   read,answ_exxf
   answ_exxf=strcompress(answ_exxf,/remove_all)
   answ_exxf=strmid(answ_exxf,0,1)
   answ_exxf=strupcase(answ_exxf)
   rrr=where(answ_exxf eq ['X','C'])
   if(rrr[0] eq (-1)) then begin
    print,' '
    print,'Please, answer ''c'' or ''x''!!!!'
    print,' '
    goto,l_eorc_ac
   endif
   if(answ_exxf eq 'X') then begin 
    print,'Exiting ...'
    goto,eocmffffs
   endif
  endif
  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'
  rrr=' '
  read,rrr,prompt=''
  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'
  rrr=' '
  read,rrr,prompt=''
  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'
  rrr=' '
  read,rrr,prompt=''
  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,' '
  if(only_disk eq 0) then begin
   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
  endif
  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=' '
  read,rrr,prompt=''
  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.)
  ;
  if(only_disk eq 0) then begin 
   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
  endif else begin
   avg_sctlght_offlimb=0
  endelse
  spflfld=(spflfld>0.)
  loadct,0
  window,1,xs=wxdim,ys=wydim,title='flatfielded raw spectra -- real relative intensities'
  tvscl,congrid(spflfld,wxdim,wydim)
  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
    tvlct,rcol_sp,gcol_sp,bcol_sp
    window,1,xs=wxdim,ys=wydim,title='flatfielded raw spectra -- real relative intensities'
    tvscl,congrid(spflfld,wxdim,wydim)
    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
   tvlct,rcol_sp,gcol_sp,bcol_sp
   window,1,xs=wxdim,ys=wydim,title='flatfielded raw spectra -- real relative intensities'
   tvscl,congrid(spflfld1,wxdim,wydim)
   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,congrid(dcspvw,3,wxdim,wydim),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=' '
  read,rrr,prompt=''
  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!dcalib_ref_prof!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!dcalib_ref_prof!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.'
     print,' '
     print,'static calibration coefficient:'
     print,calibcoeff_abs
     print,' '
     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
 ;
 ;
 ;
 ; ************* SAVING THE REDUCED DATA *****************
 loadct,0
 ;
 l_svresults:
 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'
   if(only_disk eq 0) then begin
    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]
   endif
   if(only_disk eq 1) then begin
    dxpxl_arcsec=1.
    dypxl_arcsec=1.
    xcent_pxl=0
    ycent_pxl=0
   endif    
   if(only_disk ne 1) then begin
    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 lower horizontal wire with the slit
   endif
  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=basepath,/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
  ;
  filename2dandt,flnm_sp,date_sp,time_sp
  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
     info=strarr(12)
     info[0]='spflfld_counts ........ flatfielded spectrum, intensities in counts'
     info[1]='spflfld_calibrated .... pectrum calibrated into physical units'
     info[2]='wvl ................... vector of wavelength'
     info[3]='date_sp ............... date of the spectral observations'
     info[3]='time_sp ............... time of spectral observations in UT'
     info[4]='ypos_wires_sp ......... y-pxl-positions of two horizontal wires in the spectrum'
     info[5]='calib_info ............ general info about absolute calibration (esp. physical units)'
     info[6]='imgsj ................. slit-jaw image corresponding to the spectrum'
     info[7]='xcoords_sj_arcsec, ycoords_sj_arcsec .... vectors of x- and y-coordinates in arcsec for the slit-jaw image'
     info[8]='date_sj ............... date of the slit-jaw observation'
     info[9]='time_sj ............... time of the slit-jaw observation in UT'
     info[10]='uwire_sj_pxl ......... [x,y] pxl position of upper hair in slit-jaw image'
     info[11]='lwire_sj_pxl ......... [x,y] pxl position of lower hair in slit-jaw image'
     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,info,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,' '
   if(only_disk ne 1) then begin
    flsjout=path_out+flnm_sj+'_prep.idl'
    info1=strarr(9)
    info1[0]='IMGSJ ...... slit-jaw image'
    info1[1]='DATE_SJ .... date of SJ observation'
    info1[2]='TIME_SJ .... time of SJ observation'
    info1[1]='UPHAIRX .... x position of top hair on the slit [pxl]'
    info1[2]='UPHAIRY .... y position of top hair on the slit [pxl]'
    info1[3]='DWNHAIRX .... x position of bottom hair on the slit [pxl]'
    info1[4]='DWNHAIRY .... y position of bottom hair on the slit [pxl]'
    info1[5]='CRPIX1 .... x position of the disk center [pxl]'
    info1[6]='CRPIX2 .... y position of the disk center [pxl]'
    info1[7]='CDELT1 x-size of the pixel [arcsec]'
    info1[8]='CDELT2 y-size of the pixel [arcsec]'  
    UPHAIRX=fix(round(uwire_sj_pxl[0]))
    UPHAIRY=fix(round(uwire_sj_pxl[1]))
    DWNHAIRX=fix(round(lwire_sj_pxl[0]))
    DWNHAIRY=fix(round(lwire_sj_pxl[1]))
    CRPIX1=fix(round(xcent_pxl))
    CRPIX2=fix(round(ycent_pxl))
    CDELT1=round(dxpxl_arcsec*1.e4)/1.e4
    CDELT2=round(dypxl_arcsec*1.e4)/1.e4
    save,imgsj,date_sj,time_sj,info1,UPHAIRX,UPHAIRY,DWNHAIRX,DWNHAIRY,CRPIX1,CRPIX2,CDELT1,CDELT2,file=flsjout
   endif
  endif
  if(fltype eq 'fts') then begin
   if(sw_abscal eq 'Y') then flspout=path_out+flnm_sp+'_int.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,'SCTLDSK',scattered_light_disk,'subtracted scattered light at disk'
     sxaddpar,outhdr,'CNTARLL',twocontareas[0],'left wvl border of a continuum area selected on left from sp.l.'
     sxaddpar,outhdr,'CNTARLR',twocontareas[1],'right wvl border of a continuum area selected on left from sp.l.'
     sxaddpar,outhdr,'CNTARRL',twocontareas[2],'left wvl border of a continuum area selected on right from sp.l.'
     sxaddpar,outhdr,'CNTARRR',twocontareas[3],'right wvl border of a continuum area selected on right from sp.l.'
     sxaddpar,outhdr,'SCTLLMB',avg_sctlght_offlimb,'subtracted off-limb scattered light'
     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)'
     sxaddpar,outhdr,'COMMENT','CNTARLL, CNTARLR, CNTARRL and CNTARRR are wvl borders of two sections in continuum'
     sxaddpar,outhdr,'COMMENT',' used for estimation of the avg off-limb scattered light'
    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+'_cnt.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,' '
    ;
    if(only_disk ne 1) then begin
     flsjout=path_out+flnm_sj+'_prep.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'
    endif
    print,' '
   endif
  endif
 endif
endif
;
eocmffffs:
print,' '
print,' '
print,'  ********************************'
print,'  *   Code finished its run ...  *'
print,'  ********************************'
print,' '
print,' '
end