;---  pro memu

PRO MENU, BUT_DESC, OPTION, TITLE=TITLE, XOFFSET=XOFFSET, YOFFSET=YOFFSET, $
      ROW=ROW, BUTSIZE=BUTSIZE, BASE=BASE, BUTTONS=BUTTONS, FONT=FONT, $
      frame=frame, ypad=ypad, no_release = no_release, $
      scroll=scroll, space=space, xpad=xpad, exclusive=exclusive, $
      uvalue=uvalue, x_scroll_size=x_scroll_size, $
      y_scroll_size=y_scroll_size


parent = 0
if (not keyword_set(TITLE)) then TITLE = 'MENU'
s = size(but_desc)
value_type = s(s(0) + 1)
if ((value_type ne 1) and (value_type ne 7)) then $
    message, 'values must be a string vector or 3-d byte array.`
  if (value_type eq 1) then begin
    if (s(0) ne 3) then message, 'type byte values must be 3-d'
    n_buttons = s(3)
  endif else begin
    n_buttons = n_elements(but_desc)
endelse

; sort out the keywords
if (not keyword_set(butsize)) then butsize=100
if (not keyword_set(xoffset)) then xoffset=5
if (not keyword_set(yoffset)) then yoffset=25
if (not keyword_set(row)) then begin
    row=0
    column=1
  endif else begin
    row=1
    column=0
endelse
if (not keyword_set(font)) then font = ''
if (not keyword_set(exclusive)) then exclusive=0
if (not keyword_set(nonexclusive)) then nonexclusive=0
if (keyword_set(scroll) or keyword_set(x_scroll_size) or $
      keyword_set(y_scroll_size)) then begin
     scroll = 1;
     if (not keyword_set(x_scroll_size)) then x_scroll_size=0
     if (not keyword_set(y_scroll_size)) then y_scroll_size=0
     endif else begin
     scroll=0
endelse
if (not keyword_set(frame)) then frame = 0
if (not keyword_set(space)) then space = 0
if (not keyword_set(xpad)) then xpad = 0
if (not keyword_set(ypad)) then ypad = 0
if (not keyword_set(uvalue)) then begin
    uvalue=lindgen(n_buttons)
  endif else begin
    s = size(uvalue)
    if (s(s(0) + 2) ne n_buttons) then $
      message, 'uvalue must have the same number of elements as values'
endelse

; create the base
  if (parent eq 0) then begin
    if (scroll) then $
      base = widget_base(column=column, exclusive=exclusive, $
      frame=frame, nonexclusive=nonexclusive, row=row, scroll=scroll, $
      space=space, xpad=xpad, ypad=ypad, x_scroll_size=x_scroll_size, $
      y_scroll_size=y_scroll_size, title = title) $
    else $
      base = widget_base(column=column, exclusive=exclusive, $
      frame=frame, nonexclusive=nonexclusive, row=row, $
      space=space, xpad=xpad, ypad=ypad, title = title)
  endif else begin
    if (keyword_set(title)) then begin
      theparent = widget_base(parent, /column, /frame)
      thelabel = widget_label(theparent, value = title)
    endif else theparent = parent
    if (scroll) then $
      base = widget_base(theparent, column=column, exclusive=exclusive, $
      frame=frame, nonexclusive=nonexclusive, row=row, scroll=scroll, $
      space=space, xpad=xpad, ypad=ypad, x_scroll_size=x_scroll_size, $
      y_scroll_size=y_scroll_size) $
    else $
      base = widget_base(theparent, column=column, exclusive=exclusive, $
      frame=frame, nonexclusive=nonexclusive, row=row, $
      space=space, xpad=xpad, ypad=ypad)
  endelse

; create the buttons
buttons = lindgen(n_buttons)
if (value_type eq 1) then begin
    for i = 0, n_buttons-1 do $
      buttons(i) = widget_button(base, value=but_desc(*, *, i), $
       no_release = no_release, uvalue=uvalue(i), $
       font = font, xsize=butsize)
endif else begin
    for i = 0, n_buttons-1 do $
      buttons(i) = widget_button(base, value=but_desc(i), $
       no_release = no_release, uvalue=uvalue(i), $
       font = font, xsize=butsize)
endelse

; realizacja na ekranie

widget_control, /realize, base, tlb_set_xoffset=xoffset, tlb_set_yoffset=yoffset
event=widget_event(base)
kl= fix ( where (buttons eq event.id) + 1 )
widget_control, base, /destroy
option=kl(0)
end



;---  end of pro menu   ----------------------------------------


;--- MAIN PROGRAM   --------------------------------------------


;INSTRUCTIONS
;put all data frames you want to flat-field to one directory with dark frames and flats you want to use

;ASSUMPTIONS:
;shift between DF, FF and DATA is insignificant (negligible) - unfortunettly not true
;in one set of data (approx 1-2 hours)

;number of cameras to be analysed
;0 sj
;1 ha
;2 ca
;3 d3
;4 hb
;5 ir

;-----------------

;free memory:
dark=0
ave_lines=0
flat=0
ma_dark = 0
ma_flat = 0
ma_flat_ave = 0
ma_flat_sub = 0
med_dark = 0
nol_ma_flat_ave = 0
row_ma_flat_ave = 0


!x.style=1

cam_index = strarr(6)
litt = strarr(6)

litt(0) = 'sj'
litt(1) = 'ha'
litt(2) = 'ca'
litt(3) = 'd3'
litt(4) = 'hb'
litt(5) = 'ir'

NO_CAM = 0

repeat begin

    MENU,$
    ['SJ',$
     'HA',$
     'CA',$
     'D3',$
     'HB',$
     'IR',$
     '.:: OK ::.'],$
    LETR,TITLE='SELECT CAMERAS (one click for one camera)',BUTSIZE=635,YOFFSET=50

    IF LETR LT 7 THEN CAM_INDEX(NO_CAM) = litt(letr-1)

    IF LETR NE 7 THEN no_cam = no_cam + 1

ENDREP UNTIL LETR EQ 7


;-----------------


dir=dialog_pickfile(/DIRECTORY)
sub_dir = 'resources\'
sub_dir_name = 'resources'
file_mkdir, sub_dir_name

openw,1,dir+sub_dir+'cameras.txt'
for i=0,no_cam-1 do printf,1,cam_index(i)
close,1

;--------   DARK FRAME  --------------------------------



MENU,['CALCULATE MASTER DARK','READ MASTER DARKS FILES','QUIT'],LETR,TITLE='CONTINUE?',BUTSIZE=800
;   IF (LETR EQ 1) THEN
    IF (LETR EQ 2) THEN BEGIN
       ma_dark = fltarr (1280, 512, no_cam)
       for j=0,no_cam-1 do begin
          ma_dark(*,*,j)= readfits(dir+sub_dir+'master_dark_'+cam_index(j)+'.fit')
       endfor
       GOTO,AFTER_DARK
    ENDIF
    IF (LETR EQ 3) THEN GOTO,KONIEC


;finding dark frames:

;determination of number of darks
no_D = intarr(no_cam)
filelist_D = FINDFILE(dir+'*'+cam_index(0)+'_D*.fit')
no_D(0) = size(filelist_D,/N_ELEMENTS)

;allocate array to filelist
filelist_D = strarr(no_D(0),no_cam)


;finding files
for j=0,no_cam-1 do begin
    ;check number of darks
    no_D(j) = size(FINDFILE(dir+'*'+cam_index(j)+'_D_*.fit'),/N_ELEMENTS)
    filelist_D(0:NO_D(J)-1,j) = FINDFILE(dir+'*'+cam_index(j)+'_D_*.fit')
endfor

;number of dark frames

;for j=0,no_cam-1 do begin;
    ;no_D(j) = size(filelist_D(*,j),/N_ELEMENTS)
;endfor

;control if number of darks are equal for all cameras:
control=0

for j=1,no_cam-1 do begin
    if no_D(0) ne no_D(j) then control=control+1
endfor


if control eq 0 then begin
    MENU,['YES','NO'],LETR,TITLE='NUMBER OF DARKS:'+ string(no_D(0)) + '   DO YOU WANT TO CONTINUE?',BUTSIZE=800
;   IF (LETR EQ 1) THEN
    IF (LETR EQ 2) THEN GOTO,KONIEC
    endif else begin
    MENU,['YES','NO'],LETR,TITLE='NUMBER OF DARKS ARE NOT EQUAL IN ALL CAMERAS, DO YOU WANT TO CONTINUE?',BUTSIZE=800
;   IF (LETR EQ 1) THEN
    IF (LETR EQ 2) THEN GOTO,KONIEC
    endelse

;

;reading darks:


MENU,['FIND BAD PIXELS','READ BAD PIXEL FILES','QUIT'],LETR,TITLE='CONTINUE?',BUTSIZE=800
    ;  IF (LETR EQ 1) THEN
       IF (LETR EQ 2) THEN GOTO,AFTER_DARK
       IF (LETR EQ 3) THEN GOTO,KONIEC

MENU,['MEAN VALUE','MEDIAN VALUE'],LETR_METHOD,TITLE='CHOOSE CALCULATING MASTER DARK FRAME METHOD',BUTSIZE=500


ma_dark = dblarr(1280,512,NO_CAM)

;------------------     LOOP FOR DARKS   -------------------------------

FOR J=0,NO_CAM-1 DO BEGIN
    ;master dark:

    dark = FLTARR(no_D(J),1280,512)

    for i=0,no_D(J)-1 do begin
          dark[i,*,*]=readfits(filelist_D[i,j],header)
          print,'reading darks',i,filelist_D[i,j]
    endfor


    case letr_METHOD of

    1:  begin  ;calculating master dark from mean value of pixels:

          for i=0,no_D(J)-1 do begin
              ma_dark(*,*,j) = ma_dark(*,*,j) + dark(i,*,*)
          endfor
          ma_dark(*,*,j) = ma_dark(*,*,j) / no_D(J)
          print,'calculating master dark [mean value]',j

       end


    2:  begin  ;calculating master dark from median value for each pixel:

          print,'calculating master dark [median value]',j
          for m=0,1279 do begin
              for n=0,511 do begin
                 ma_dark[m,n,J] = median(dark[*,m,n])
              endfor
          endfor

       end

    endcase


       window,j,ypos=j*100,xpos=j*50,xsize=640,ysize=256,title=cam_index(j)+'master dark'
       tvscl,rebin(ma_dark(*,*,j),640,256)

    ;---  detection of bad pixels

    ;caclulating median value of master darks:

    ;med_dark = dblarr(no_cam)


       med_dark = median(ma_dark(*,*,j))



;absolute devation:
;dev=dblarr(no_cam)


    dev = total(abs(ma_dark(*,*,j) - med_dark)) / (1280.*512.)




    ;cut = dblarr(no_cam)
    ;bad_pxl = intarr(no_cam,1024)

    x_ = indgen(max(ma_dark(*,*,j))-min(ma_dark(*,*,j))+1)+min(ma_dark(*,*,j))

    wait,0.5

    ;showing histograms
    window,0,ypos=0,xpos=0,xsize=640,ysize=256,title=cam_index(j)+' master dark histogram'
    plot,x_,histogram(ma_dark(*,*,j),binsize=1),psym=10,/ylog,yrange=[0.15,1000000]

    window,1,ypos=0,xpos=645,xsize=640,ysize=256,title=cam_index(j)+' master dark'
    tvscl,rebin(ma_dark(*,*,j),640,256)


    ;-------   INTERACTIVE MODULE TO CHOOSE BAD PIXEL CUTTING LEVEL

    loadct,3

    ;FOR SJ DARK
    REPEAT BEGIN
       MENU,$
       [' 3 SIGMA='+STRING(DEV*3.)+' ADU',$
        ' 5 SIGMA='+STRING(DEV*5.)+' ADU',$
        ' 7 SIGMA='+STRING(DEV*7.)+' ADU',$
        ' 9 SIGMA='+STRING(DEV*9.)+' ADU',$
        '11 SIGMA='+STRING(DEV*11.)+' ADU',$
        '13 SIGMA='+STRING(DEV*13.)+' ADU',$
        '15 SIGMA='+STRING(DEV*15.)+' ADU',$
        '17 SIGMA='+STRING(DEV*17.)+' ADU',$
        'OK'],$
       LETR,TITLE=cam_index(j) +' ABSOLUTE DEVIATION:'+STRING(DEV)+'  CHOOSE CUTTING LEVEL',BUTSIZE=635,YOFFSET=275

       CASE LETR OF
         1:  cut = dev*3.
         2:  cut = dev*5.
         3:  cut = dev*7.
         4:  cut = dev*9.
         5:  cut = dev*11.
         6:  cut = dev*13.
         7:  cut = dev*15.
         8:  cut = dev*17.
         9:
       ENDCASE

       ;subscripts(one-dimentional) of bad pixels on master dark frame:
       bad_pxl=where(ma_dark(*,*,j) GT med_dark+cut OR ma_dark(*,*,j) LT med_dark-cut, no_bad_pxl)


       window,0,ypos=0,xpos=0,xsize=640,ysize=256,title=cam_index(j)+' master dark histogram, no. of bad pixels:'+string(no_bad_pxl)
       plot,x_,histogram(ma_dark(*,*,j),binsize=1),psym=10,/ylog,yrange=[0.15,1000000]
       PLOTS,[med_dark-cut,med_dark-cut],[0.01,1000000],/Data,color=150
       PLOTS,[med_dark+cut,med_dark+cut],[0.01,1000000],/Data,color=150
    ENDREP UNTIL LETR EQ 9

    ;trick to format output data:
    ;to write to file in column order
    bad_write=lonarr(1,size(bad_pxl,/n_elements))
    bad_write(0,*)=bad_pxl

    openw,1,dir+sub_dir+'bad_pxl_'+cam_index(j)+'.txt'
    printf,1,bad_write
    close,1

    ;---   replacing bad pixels value with mean of 4 nearby pixels wtih indices: [-1,-1],[-1,+1],[+1,-1],[+1,+1]

    ;bad_col = bad_pxl mod 1280
    ;bad_row = bad_pxl / 1280
    ;;bad_sj = [[bad_col_sj],[bad_row_sj]]

    ;ma_dark[bad_col,bad_row,j] = $
    ;(ma_dark(bad_col-1,bad_row+1,j)+$
    ;ma_dark(bad_col-1,bad_row-1,j)+$
    ;ma_dark(bad_col+1,bad_row+1,j)+$
    ;ma_dark(bad_col+1,bad_row-1,j))/4.



;for j=0,no_cam-1 do begin
;   window,j,ypos=j*100,xpos=j*50,xsize=640,ysize=256,title=cam_index(j)+'master dark'
;   tvscl,rebin(ma_dark(*,*,j),640,256)
;endfor

;------  END OF INTERACTIVE MODULE



;for j=0,no_cam-1 do begin
    writefits,dir+sub_dir+'master_dark_'+cam_index(j)+'.fit',ma_dark(*,*,j)
;Endfor


endfor

;-------------    END OF LOOP FOR DARKS   -------------------------

AFTER_DARK:

;----------   END OF DARK FRAME MODULE   -------------------


;free memory:
dark=0
med_dark = 0
ma_flat_ave = fltarr(1280, 512, no_cam)


;----------   FLAT FIELD MODULE   --------------------------


MENU,['CALCULATE MASTER FLAT','READ MASTER FLAT FILES','QUIT'],LETR,TITLE='CONTINUE?',BUTSIZE=800
;   IF (LETR EQ 1) THEN
    IF (LETR EQ 2) THEN BEGIN
       for j=0,no_cam-1 do begin
         ma_flat_ave(*,*,j)= readfits(dir+sub_dir+'master_flat_'+cam_index(j)+'.fit',header)
       endfor
       GOTO,AFTER_FLAT
    ENDIF
    IF (LETR EQ 3) THEN GOTO,KONIEC

;find flat field frames:


;number of flat frames
no_F = intarr(6)

for j=0,no_cam-1 do begin
    filelist_F = FINDFILE(dir+'*'+cam_index(j)+'_F*.fit')
    no_F(j) = size(filelist_F,/N_ELEMENTS)
endfor

;control if number of flats are equal for all cameras:
control=0

for j=1,no_cam-1 do begin
    if no_F(0) ne no_F(j) then control=control+1
endfor

filelist_F = strarr(6,max(no_F))

if control ne 0 then begin
    for i=0,no_cam-1 do begin
       print,'number of FF(',cam_index(i),')',no_F(i)
    endfor
    print,'Delete group of FF which are not complete by all cameras'
    ;stop
endif

for j=0,no_cam-1 do begin
    temp_filelist = FINDFILE(dir+'*'+cam_index(j)+'_F*.fit')
    temp_size = N_ELEMENTS(temp_filelist)
    FOR I=0,temp_SIZE-1 DO BEGIN
       filelist_F(j,I) = temp_filelist(I)
    ENDFOR
endfor




if control eq 0 then begin
    MENU,['YES','NO'],LETR,TITLE='NUMBER OF FLATS:'+ string(no_F(0)) + '   DO YOU WANT TO CONTINUE?',BUTSIZE=600
;   IF (LETR EQ 1) THEN
    IF (LETR EQ 2) THEN GOTO,KONIEC
    endif else begin
    MENU,['YES','NO'],LETR,TITLE='NUMBER OF FLATS ARE NOT EQUAL IN ALL CAMERAS, DO YOU WANT TO CONTINUE?',BUTSIZE=600
;   IF (LETR EQ 1) THEN
    IF (LETR EQ 2) THEN GOTO,KONIEC
    endelse










;master flat:
ma_flat = fltarr(1280,512,no_cam)
MENU,['MEAN VALUE','MEDIAN VALUE'],LETR,TITLE='CHOOSE CALCULATING MASTER FLAT FIELD',BUTSIZE=500

;reading flats:
; camera by camera:

for j=0,no_cam-1 do begin

    flat = intarr(no_F(j),1280,512)

    for i=0,no_F(j)-1 do begin
       IF filelist_F[j,i] NE '' THEN BEGIN
         flat[i,*,*]=readfits(filelist_F[j,i],header)
       ENDIF
       if i mod 10 eq 0 then begin
         print,'reading flats',j
       endif
    endfor
    print,'complete'

    case letr of

    1:  begin  ;calculating master flat from mean value of pixels:

          for i=0,no_F(0)-1 do begin
              ma_flat(*,*,j) = ma_flat(*,*,j) + flat[i,*,*]
          endfor
          ma_flat(*,*,j) = ma_flat(*,*,j) / no_F(j)
          print,'calculating master flat [mean value]',j

       end

    2:  begin  ;calculating master flat from median value for each pixel:

         print,'calculating master flat [median value]',j
          for m=0,1279 do begin
              for n=0,511 do begin
                 ma_flat[m,n,j] = median(flat[*,m,n])
              endfor
          endfor
       end

    endcase

endfor














;---


;showing master flats:




for j=0,no_cam-1 do begin
    window,j,ypos=j*100,xpos=j*50,xsize=640,ysize=256,title=cam_index(j)+' master flat'
    tvscl,rebin(ma_flat(*,*,j),640,256)
endfor

MENU,['YES','QUIT'],LETR,TITLE='CONTINUE?',BUTSIZE=800
;   IF (LETR EQ 1) THEN
    IF (LETR EQ 2) THEN GOTO,KONIEC

FOR I=0,NO_CAM-1 DO WDELETE,I


;subtracting darks from flats:
ma_flat_sub=dblarr(1280,512,no_cam)

for j=0,no_cam-1 do begin
    ma_flat_sub(*,*,j) = ma_flat(*,*,j) - ma_dark(*,*,j)
endfor

;normalization of flats:
ma_flat_ave = dblarr(1280,512,no_cam)

for j=0,no_cam-1 do begin
    ma_flat_ave(*,*,j) = ma_flat_sub(*,*,j)/ median(ma_flat_sub(*,*,j))
endfor


wait,1.0




;---   replacing bad pixels value with mean of 4 nearby pixels with indices: [-1,-1],[-1,+1],[+1,-1],[+1,+1]

for j=0,no_cam-1 do begin

    ;another trick
    ;in this case first calculate size of file:
    CLOSE,/ALL
    ON_IOERROR, B222
    A1=' '
    OPENR,1,dir+sub_dir+'bad_pxl_'+cam_index(j)+'.txt'
    FOR I=0,999 DO READF,1,A1
    B222: no_bad=I
    CLOSE,1
    ;PRINT,no_bad
    ON_IOERROR,NULL

    ;then creating array with corresponding size:
    bad_pxl = lonarr(no_bad)

    ;finally read the file again:
    openr,1,dir+sub_dir+'bad_pxl_'+cam_index(j)+'.txt'
    readf,1,bad_pxl
    close,1


    if size(bad_pxl,/n_elements) ne 1 and bad_pxl(0) ne -1 then begin

       bad_col = bad_pxl mod 1280
       bad_row = bad_pxl / 1280
       ;bad_sj = [[bad_col_sj],[bad_row_sj]]
       print,'no_bad:',no_bad,'   cam_index:',j
       ;print,'bad_col',bad_col
       ;print,'bad_row',bad_row
       jj=replicate(j,no_bad)
       print,'median values of FF:',median(ma_flat(*,*,j))
       print,'values of bad pxls on FF:',ma_flat_ave(bad_col,bad_row,jj)
       print,'values of bad pxls on DF:',ma_dark(bad_col,bad_row,jj)

        ma_flat_ave[bad_col,  bad_row,  jj] = $
       (ma_flat_ave(bad_col-1,bad_row+1,jj) +$
        ma_flat_ave(bad_col-1,bad_row-1,jj) +$
        ma_flat_ave(bad_col+1,bad_row+1,jj) +$
        ma_flat_ave(bad_col+1,bad_row-1,jj))/4.

    endif else begin
       print,'   cam_index:',j
       print,'no bad pixels detected'
    endelse


endfor

;tu sie zatrzymalem
;postprawdzac usrednianie bad pxl


for j=0,no_cam-1 do begin  ;writing normalized, bad pixel removed master flats:
    writefits,dir+sub_dir+'master_flat_'+cam_index(j)+'.fit',ma_flat_ave(*,*,j)
endfor


AFTER_FLAT:

;showing normalized master flats

for j=0,no_cam-1 do begin
    window,j,ypos=j*100,xpos=j*50,xsize=640*2,ysize=256*2,title=cam_index(j)+' normalized master flat'
    tvscl,ma_flat_ave(*,*,j)
endfor

;----------   END OF FLAT FIELD MODULE



MENU,['YES','QUIT'],LETR,TITLE='CONTINUE?',BUTSIZE=800
;   IF (LETR EQ 1) THEN
    IF (LETR EQ 2) THEN GOTO,KONIEC

for i=0,32 do wdelete,i


loadct,3

;--------    SPECTRAL LINE CORRECTION MODULE  ----------------


line_pos_coef = dblarr(3,no_cam)
line_pos_poly = dblarr(512,no_cam)
x_coord = findgen(512)
pxl_shift = dblarr(no_cam) ; variable for value of max shift for each image

;---  determination of curvature and/or rotation of spectra (and slit jaw image) on ccd images
;click on image to mark line region

for i=0,no_cam-1 do begin

SECOND_LINE = 0

    REPEAT BEGIN

       window,6,ypos=0,xpos=0,xsize=1280,ysize=512,title='mark slit/line position on the image'
       tvscl,ma_flat_ave(*,*,i)

       CURSOR, x, Y, /device, /down
       x_line_pos_1 = x;*2
       PLOTS,[x,x],[0,512],/device,color=255

       CURSOR, x, Y, /device, /down
       x_line_pos_2 = x;*2
       PLOTS,[x,x],[0,512],/device,color=255

       wait,0.4

       x_line_pos_left=x_line_pos_1<x_line_pos_2
       x_line_pos_right=x_line_pos_2>x_line_pos_1

       line_pos = dblarr(512)

       for j=0,511 do begin
         line = where( ma_flat_ave[x_line_pos_left:x_line_pos_right,j,i] eq  $
                  min( ma_flat_ave[x_line_pos_left:x_line_pos_right,j,i] ) ) $
                  + x_line_pos_left
         line_pos(j) = line(0)
       endfor

       line_pos_1 = dblarr(512)
       poly_error = dblarr(512)
       y_coord    = indgen(1280)
       chisq_tab  = fltarr(512)
       yerror_tab = fltarr(512)

       for j=0,511 do begin
         f=poly_fit(y_coord[line_pos(j)-3:line_pos(j)+3],ma_flat_ave(line_pos(j)-3:line_pos(j)+3,j,i),2,/double,sigma=sigma,chisq=chisq, yerror=yerror)
         chisq_tab(j)  = chisq
         yerror_tab(j) = yerror
         line_pos_1(j) = -f(1)/(2*f(2))
         poly_error(j) = sqrt(  (sigma(2) * f(1) / ( 2.*f(2)^2. ))^2.  +  (sigma(1) / ( 2.*f(2) ))^2. )
       endfor

       line_pos_1 = line_pos_1 > x_line_pos_left
       line_pos_1 = line_pos_1 < x_line_pos_right

       ;quadratic approximation:
       line_pos_coef(*,i) = poly_fit(x_coord,line_pos_1,2,/double,yfit=y_pos_line)  ;,MEASURE_ERRORS=poly_error,;;;;;; <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
       line_pos_poly(*,i) = line_pos_coef(0,i) +  x_coord*line_pos_coef(1,i) + x_coord*x_coord*line_pos_coef(2,i)

       print,'CAMERA:  ',i
       print,'COEFF:   ',line_pos_coef(*,i)

       window,0,ypos=560,xpos=650,xsize=640,ysize=256,title='slit/line position along vertical axis'
       plot,x_coord,line_pos_1,psym=6,symsize=0.3,/ynozero,ystyle=2,yrange=[x_line_pos_left,x_line_pos_right]
       oplot,x_coord,line_pos_poly(*,i),color=150,thick=3
       pxl_shift(i) = 512.*line_pos_coef(1,i)+512.*512.*line_pos_coef(2,i)
       print,'pxl_shift       (top -- bottom diff): ', line_pos_poly(511,i) - line_pos_poly(0,i)
       print,'pxl_shift (maximum difference): ', max(line_pos_poly(*,i)) - min(line_pos_poly(*,i))

       MENU,$
       ['AGAIN',$
        'SECOND LINE', $
        'OK'],$
       LETR,TITLE='POLY FIT APPROXIMATION',BUTSIZE=635,YOFFSET=560

       IF LETR EQ 2 THEN BEGIN
         IF SECOND_LINE EQ 1 THEN BEGIN
          MENU,['OK'],RRR,TITLE='OH, NO MY DARLING. YOU COULD NOT SELECT MORE THAN TWO LINES. SORRY.',BUTSIZE=700,YOFFSET=360
          LETR = 3
         ENDIF ELSE BEGIN
          SECOND_LINE = 1
          FIRST_LINE_POS_POLY = line_pos_poly(*,i)
         ENDELSE
       ENDIF

    ENDREP UNTIL LETR EQ 3

    IF SECOND_LINE EQ 1 THEN BEGIN
       line_pos_poly(*,i) = (FIRST_LINE_POS_POLY + line_pos_poly(*,i)) / 2.
    ENDIF

endfor



;---  shifting individual rows of pixels in horizontal direction to make spectral lines strictly vertical:
;note that from this moment pixels near the edge of images are not correct because they came from shifting procedure

row_calc = dblarr(1280*100.)
row_ma_flat_ave = dblarr(1280,512,no_cam)

for i=0,no_cam-1 do begin
    print,'shifting rows...'
       for j=0,511 do begin
       row_calc=rebin(ma_flat_ave(*,j,i),1280*100.)
          ;j_float = double(j)
          ;row_shift=shift    ( row_calc, round(-100.*   (j_float*line_pos_coef(1,i)  )  + j_float*j_float*line_pos_coef(2,i)   )  );<<<<<<<<<<<<<<<<<<
       row_shift=shift   ( row_calc, round(-100.* (line_pos_poly(j,i)-line_pos_poly(0,i) ) )  )
       row_write=rebin(row_shift,1280.)
       row_ma_flat_ave(*,j,i)=row_write
    endfor
    print,'camera',i,'   complete'
endfor


;--------    END OF SPECTRAL LINE CORRECTION MODULE   ---------------

;---------   SPECTRAL LINE ELIMINATION MODULE   ----------------

;calculating mean/median value of each column

ave_col=dblarr(1280,no_cam)

for i=1,no_cam-1 do begin
    for j=0,1279 do begin
       ave_col(j,i) = median(row_ma_flat_ave(j,*,i))
    endfor
endfor

;creating 1280x512 pxl matrix with mean/median value of each column replicated to whole column

ave_lines = dblarr(1280,512,no_cam)

for i=1,no_cam-1 do begin
    for  j=0,1279 do begin
       ave_lines(j,*,i) = replicate(ave_col(j,i),1,512)
    endfor
endfor

;------

for j=1,no_cam-1 do begin
    window,j,ypos=j*100,xpos=j*50,xsize=1280/2,ysize=512/2,title=cam_index(j)+' ave lines matrices'
    tvscl,rebin(ave_lines(*,*,j),1280/2,512/2)
endfor


;--------


;elimination lines from spectral master flats

nol_ma_flat_ave = dblarr(1280,512,no_cam)

for i=1,no_cam-1 do begin
    nol_ma_flat_ave(*,*,i) = row_ma_flat_ave(*,*,i) / ave_lines(*,*,i)
endfor

;only for comfort,no_cam = 0 means sj-image
nol_ma_flat_ave(*,*,0)=row_ma_flat_ave(*,*,0)

for j=0,no_cam-1 do begin
    window,j,ypos=j*100,xpos=j*50,xsize=1280,ysize=512,title=cam_index(j)+' normalized master flat after line elimination'
    tvscl,nol_ma_flat_ave(*,*,j)
endfor

MENU,['YES','QUIT'],LETR,TITLE='CONTINUE?',BUTSIZE=800
;   IF (LETR EQ 1) THEN
    IF (LETR EQ 2) THEN GOTO,KONIEC

FOR I=0,NO_CAM-1 DO WDELETE,I

for j=0,no_cam-1 do begin
    writefits,dir+sub_dir+'nol_flat_'+cam_index(j)+'.fit',nol_ma_flat_ave(*,*,j)
endfor


;;--------  HAIRS POSITION CORRECTION MODULE  ------------------------

y_coord = findgen(1280)

hair_pos_up_coef = dblarr(3,6)
hair_pos_up_poly = dblarr(1280,6)

hair_pos_down_coef = dblarr(3,6)
hair_pos_down_poly = dblarr(1280,6)

for i=0,no_cam-1 do begin

    REPEAT BEGIN

       window,6,ypos=0,xpos=0,xsize=1280,ysize=512,title='mark hairs position on the image, first top then bottom'
       tvscl,nol_ma_flat_ave(*,*,i)

       CURSOR, x, Y, /device, /down
       y_hair_pos_1 = y;*2
       print,y
       PLOTS,[0,1279],[y,y],/device,color=255

       CURSOR, x, Y, /device, /down
       y_hair_pos_2 = y;*2
       print,y
       PLOTS,[0,1279],[y,y],/device,color=255

       y_hair_1_pos_down=y_hair_pos_1<y_hair_pos_2
       y_hair_1_pos_up=y_hair_pos_2>y_hair_pos_1


       CURSOR, x, Y, /device, /down
       y_hair_pos_1 = y;*2
       print,y
       PLOTS,[0,1279],[y,y],/device,color=255

       CURSOR, x, Y, /device, /down
       y_hair_pos_2 = y;*2
       print,y
       PLOTS,[0,1279],[y,y],/device,color=255

       y_hair_2_pos_down=y_hair_pos_1<y_hair_pos_2
       y_hair_2_pos_up=y_hair_pos_2>y_hair_pos_1

wait,0.4

       hair_pos_up = dblarr(1280)
       hair_pos_down = dblarr(1280)

       for j=0,1279 do begin
         hair = where( nol_ma_flat_ave[j,y_hair_1_pos_down:y_hair_1_pos_up,i] eq  $
                  min( nol_ma_flat_ave[j,y_hair_1_pos_down:y_hair_1_pos_up,i] ) ) $
                  + y_hair_1_pos_down
         ;if size(hair,/N_elements) ne 1 then hair_tot = total(hair) / size(hair,/N_elements) else hair_tot = hair
         hair_pos_up(j) = hair(0)
       endfor

       for j=0,1279 do begin
         hair = where( nol_ma_flat_ave[j,y_hair_2_pos_down:y_hair_2_pos_up,i] eq $
                min( nol_ma_flat_ave[j,y_hair_2_pos_down:y_hair_2_pos_up,i] ) ) $
                + y_hair_2_pos_down
         ;if size(hair,/N_elements) ne 1 then hair_tot = total(hair) / size(hair,/N_elements) else hair_tot = hair
         hair_pos_down(j) = hair(0)
       endfor

       hair_pos_up_1 = dblarr(1280)
       hair_pos_down_1 = dblarr(1280)
       poly_error = dblarr(1280)
       y_coord    = findgen(1280)
       chisq_tab  = fltarr(1280)
       yerror_tab = fltarr(1280)

       for j=0,1279 do begin
         f=poly_fit(y_coord[hair_pos_up(j)-2:hair_pos_up(j)+2],ma_flat_ave(j, hair_pos_up(j)-2:hair_pos_up(j)+2,i),2,/double,sigma=sigma,chisq=chisq, yerror=yerror)
         chisq_tab(j)  = chisq
         yerror_tab(j) = yerror
         hair_pos_up_1(j) = -f(1)/(2*f(2))
         poly_error(j) = sqrt(  (sigma(2) * f(1) / ( 2.*f(2)^2. ))^2.  +  (sigma(1) / ( 2.*f(2) ))^2. )
       endfor

       hair_pos_up_1 = hair_pos_up_1 > y_hair_1_pos_down
       hair_pos_up_1 = hair_pos_up_1 < y_hair_1_pos_up



       for j=0,1279 do begin
         f=poly_fit(y_coord[hair_pos_down(j)-2:hair_pos_down(j)+2],ma_flat_ave(j, hair_pos_down(j)-2:hair_pos_down(j)+2,i),2,/double,sigma=sigma,chisq=chisq, yerror=yerror)
         chisq_tab(j)  = chisq
         yerror_tab(j) = yerror
         hair_pos_down_1(j) = -f(1)/(2*f(2))
         poly_error(j) = sqrt(  (sigma(2) * f(1) / ( 2.*f(2)^2. ))^2.  +  (sigma(1) / ( 2.*f(2) ))^2. )
       endfor

       hair_pos_down_1 = hair_pos_down_1 > y_hair_2_pos_down
       hair_pos_down_1 = hair_pos_down_1 < y_hair_2_pos_up


       ;quadratic approximation, with included influence of shifted pixels on edges (these pixels are excluded to approximation)
       hair_pos_up_coef(*,i) =   poly_fit(y_coord[abs(pxl_shift(i)):1279-abs(pxl_shift(i))],hair_pos_up_1[abs(pxl_shift(i)):1279-abs(pxl_shift(i))],2,/double)
       hair_pos_up_poly(*,i) =   hair_pos_up_coef(0,i)   +  y_coord*hair_pos_up_coef(1,i)   + y_coord*y_coord*hair_pos_up_coef(2,i)

       hair_pos_down_coef(*,i) = poly_fit(y_coord[abs(pxl_shift(i)):1279-abs(pxl_shift(i))],hair_pos_down_1[abs(pxl_shift(i)):1279-abs(pxl_shift(i))],2,/double)
       hair_pos_down_poly(*,i) = hair_pos_down_coef(0,i) +  y_coord*hair_pos_down_coef(1,i) + y_coord*y_coord*hair_pos_down_coef(2,i)

       print,'CAMERA:  ',i
       print,'UP        HAIR COEFF:   ',hair_pos_up_coef(*,i)
       print,'DOWN HAIR COEFF:   ',hair_pos_down_coef(*,i)

       window,0,ypos=590,xpos=0,xsize=640,ysize=256,title='top hair position along horizontal axis'
       plot,y_coord,hair_pos_up_1,psym=6,symsize=0.3,/ynozero,ystyle=2,yrange=[y_hair_1_pos_down,y_hair_1_pos_up] ;here _1_ means upper hair ;)
       oplot,y_coord,hair_pos_up_poly(*,i),color=150,thick=3

       window,1,ypos=590,xpos=645,xsize=640,ysize=256,title='bottom hair position along horizontal axis'
       plot,y_coord,hair_pos_down_1,psym=6,symsize=0.3,/ynozero,ystyle=2,yrange=[y_hair_2_pos_down,y_hair_2_pos_up] ;here _2_ means down hair ;)
       oplot,y_coord,hair_pos_down_poly(*,i),color=150,thick=3



       MENU,$
       ['AGAIN',$
        'OK', $
        'XLOADCT AND AGAIN'],$
       LETR,TITLE='POLY FIT APPROXIMATION',BUTSIZE=635,YOFFSET=560

       IF LETR EQ 3 THEN XLOADCT, BLOCK=1

    ENDREP UNTIL LETR EQ 2

endfor



wait,0.5

for i=0,32 do wdelete,i


hair_pos_poly  = (hair_pos_up_poly + hair_pos_down_poly) / 2.

;---  shifting individual column of pixels in vertical direction to make hairs strictly horizontal:
;note that from this moment pixels near the edge of images are not correct because they came from shifting procedure

;--------
    col_calc = dblarr(512*20.)
    col_nol_ma_flat_ave = dblarr(1280,512,no_cam)
    ;hair_pos_coef = (hair_pos_up_coef + hair_pos_down_coef) / 2. ; make mean coefficient from 2 hairs' positions
    pxl_shift_hor = dblarr(no_cam)

    for i=0,no_cam-1 do begin
       ff=0

       for j=0,1279 do begin
         j_float = double(j)

                   ;to make 1-D array
         col_calc=rebin(total(nol_ma_flat_ave(j,*,i),1),512*20)

         scale_shift = round(20.*(  hair_pos_poly(0,i) -  hair_pos_poly(j,i) ) )
         ;scale_shift = round(    (                              - (                         +                                  ) ) )
         pxl_shift_hor(i) = scale_shift / 20.

         ff=ff+1
         if ff eq 100 then begin
          ;print,'scale:',scale,'  hair_width',hair_width
          print,'j_float:',j_float,'  scale_shift',scale_shift, '  pxl_shift_hor(i):', pxl_shift_hor(i)
          ff=0
         endif

         col_shift=shift(col_calc,scale_shift)


         col_cut = dblarr(512.*20.)
         col_cut(*) = 1. ;;  to prevent dividing by 0

         ;if scale le 1. then col_cut(0:size(col_shift,/n_elements)-1) = col_shift
         ;if scale gt 1. then col_cut = col_shift(0:512*20-1)
         col_cut = col_shift(0:512*20.-1)

         ;print,20.*j*line_pos_coef(1,i)
         col_write=rebin(col_cut,512.)
         col_nol_ma_flat_ave(j,*,i)=col_write
       endfor
       ;stop
       print,'camera',i,'   complete'
    endfor

;--------


wait,0.5


;-------- END OF HAIRS POSITION CORRECTION MODULE  ------------------


;showing normalized, bad pixel removed, corrected line and hairs positions master flats:

for j=0,no_cam-1 do begin
    window,j,ypos=j*100,xpos=j*50,xsize=640*2,ysize=256*2,title=cam_index(j)+' normalized master flat'
    tvscl,col_nol_ma_flat_ave(*,*,j)
endfor

MENU,['YES','QUIT'],LETR,TITLE='CONTINUE?',BUTSIZE=800
;   IF (LETR EQ 1) THEN
    IF (LETR EQ 2) THEN GOTO,KONIEC

FOR I=0,NO_CAM-1 DO WDELETE,I





;---   CALCULATING MEDIAN VALUE OF EACH ROW TO PRECISE FF DUE TO EXISTENCE OF HORIZONTAL FEATURES  ----

ave_row=dblarr(512,no_cam)
ave_row(*,*) = 1. ; to safety, if some element equal 0 then is dividing by 0, so it is repraced by 1 value

for i=1,no_cam-1 do begin ;from i=1 means without sj image
    for j=0,511 do begin
       ave_row(j,i) = median(col_nol_ma_flat_ave(*,j,i))
    endfor
endfor

;creating an array of mean value of each row replicared to whole row:

ave_row_array = dblarr(1280,512,no_cam)
;ave_row_array(*,*,*) = 1.

for i=1,no_cam-1 do begin
    for  j=0,511 do begin
       ave_row_array(*,j,i) = replicate(ave_row(j,i),1,1280)
       if ave_row(j,i) eq 0 then ave_row_array(*,j,i) = 1.  ;;  replicate by artificial value eq 1. to prevent dividing by 0
    endfor
endfor

for j=1,no_cam-1 do begin
    window,j,ypos=j*100,xpos=j*50,xsize=1280,ysize=512,title=cam_index(j)+' slit flat'
    tvscl,ave_row_array(*,*,j)
endfor

MENU,['YES','QUIT'],LETR,TITLE='CONTINUE?',BUTSIZE=800
;   IF (LETR EQ 1) THEN
    IF (LETR EQ 2) THEN GOTO,KONIEC

FOR I=0,NO_CAM-1 DO WDELETE,I

for j=1,no_cam-1 do begin
    writefits,dir+sub_dir+'row_flat_'+cam_index(j)+'.fit',ave_row_array(*,*,j)
endfor



;---   END OF CALCULATING MEDIAN VALUE   --------------------------

;dividing nol flat by row flat
;now nol flat not consists row flat  !!!
hard_flat = dblarr(1280,512,no_cam)

for j=1,no_cam-1 do begin
    hard_flat(*,*,j) = col_nol_ma_flat_ave(*,*,j) / ave_row_array(*,*,j)
endfor

;for comfort:
hard_flat(*,*,0) = col_nol_ma_flat_ave(*,*,0)

for j=1,no_cam-1 do begin
    window,j,ypos=j*100,xpos=j*50,xsize=1280,ysize=512,title=cam_index(j)+' hard flat'
    tvscl,hard_flat(abs(pxl_shift(j))+1:1280-abs(pxl_shift(j)-1),abs(pxl_shift_hor(j)):512-abs(pxl_shift_hor(j))-1,j) ;displaying
endfor

for j=0,no_cam-1 do begin
    writefits,dir+sub_dir+'hard_flat_'+cam_index(j)+'.fit',hard_flat(*,*,j)
endfor


;writing coefficients of lines and hairs positions:

    openw,1,dir+sub_dir+'coefficients_hair.txt'
;   for i=0,no_cam-1 do printf,1,hair_pos_coef(*,i),format='(2(F15.6))'
    for i=0,no_cam-1 do printf,1,hair_pos_up_coef(*,i),format='(3(F20.12))'
    for i=0,no_cam-1 do printf,1,hair_pos_down_coef(*,i),format='(3(F20.12))'
    close,1

    openw,1,dir+sub_dir+'coefficients_line.txt'
    for i=0,no_cam-1 do printf,1,line_pos_coef(*,i),format='(3(F20.12))'
    close,1

;---------   END OF SPECTRAL LINE ELIMINATION MODULE   ----------------


print,'SUCCES'

MENU,['OK'],LETR_METHOD,TITLE='The end of job of --- 0_FLAT --- program.   Execute --- 1_DATA ---  program',BUTSIZE=600, YOFFSET=400, XOFFSET=200


KONIEC:

print,'END OF JOB'

end