PRO DIG_MAIN4,progver,mereni,sphe,sphg,sphcor,spgcor,doppl, $
              sphe1,sphe2,sph1cor,sph2cor, $
              coel_stat,ydiam,ydiamc,runmo=RUNMO
;
; Version 5 - correction for horizontal strips included
;         H-alpha Dopplergram calculation
;         Spectroheliograms in H-alpha wings -0.5A, +0.5A
;
; Executes the processing of Coimbra spectroheliograph data. Called by ADRIDL_4.
; This version is written for Windows XP in IDL v.6.2
;
; input:  progver - program version name (string)
;         mereni - path to data diretory
; output:
;         sphe - raw spectroheliogram (centre line)
;         sphg - raw spectroheliogram (continuum)
;         sphcor - rectified spectroheliogram, no rotation (centre line)
;         spgcor - rectified spectroheliogram, no rotation (continuum)
;       the next is only for the case of H-alpha, otherwise set to 0:
;         doppl - Dopplergram, no rotation
;         sphe1 - raw spectroheliogram (H-alpha centre -0.5A)
;         sphe2 - raw spectroheliogram (H-alpha centre +0.5A)
;         sph1cor - rectified spectroheliogram, no rotation (-0.5A)
;         sph2cor - rectified spectroheliogram, no rotation (+0.5A)
; internal output:
;         coel_stat - coelostat position (string)
;         ydiam,ydiamc - solar disk diameter (line, continuum) in pixels
; keyword: runmo - run mode: 1 = single, 2 = multiple,
;                            3 = repeat preceding single run with corrections
; Note: rotated spectroheliograms and Dopplergram are NOT program
;        parameters but they are passed directly to output routines
;
; Copyright M. Klvana and M. Sobotka, ASU AVCR Ondrejov, Czech Rep., 2014
; Updated to version 5 by M. Sobotka, 2017
;

;************ SETTING OF OUTPUT PATH **********************************************

pathout='C:\12B_PIC_II\'  ;parent output data directory

;**********************************************************************************

; initial setting for strips correction mode
strip='YES'

;DIG_CORRECTIONS2 calling
if runmo eq 3 then DIG_CORRECTIONS2,mereni,coel_stat,coel,ut_cor,strip

; getting date and time
if runmo eq 3 then  $
  DIG_START,mereni,nsp,xpix,ypix,rok,mesic,den,tfirst,tlast,ut_cor=ut_cor else $
  DIG_START,mereni,nsp,xpix,ypix,rok,mesic,den,tfirst,tlast

DIG_DISK,mereni,middlesp,noframes,speed

L0B0P0,rok,mesic,den,(tfirst+tlast)/7200.,L0,B0,P0,R0

; rotation angle calculation
if runmo eq 3 then  $
  SPHE_ROT_ANGLE,rok,mesic,den,tfirst,P0,drot,coel_stat,coel=coel else  $
  SPHE_ROT_ANGLE,rok,mesic,den,tfirst,P0,drot,coel_stat

DIG_FLAT2,mereni,xpix,ypix,dc,flat,a,b,c,x0,line_name      ; flatfield

print,'Position of line:',x0
print,'Line identification: ',line_name
print,'Scanning speed: ',speed

if line_name eq 'K3' then begin           ;offsets for continuum
    offs=9.7
    line_nam2='K1'
endif else begin
    offs=16.               ;from 2017, line centre -4.0 A
    line_nam2='CN'
endelse

if runmo eq 3 then goto,skok
;-------------------------------------------------------------------------

DIG_SELECT,mereni,noframes,sellist   ;selection of spectra

; processing of selected spectra

cube=intarr(xpix,ypix,noframes)     ;volume (lambda,y,spectrum)

for i=0,noframes-1 do begin         ;loop over frames

    DIG_2,sellist(i),sp               ;reading spectra
    sp=sp-dc                          ;dark frame subtraction
    sp=sp>0
    LINERECT,sp,a,b,c,rsp             ;line curvature rectification
    rsp=rsp/flat                      ;float, flatfielded spectrum
    ssp=rsp*5                         ;re-scaling
    cube(0,0,noframes-1-i)=fix(ssp)    ;integer*2

endfor

; spectroheliograms calculation

sphe=transpose(reform(cube(round(x0),*,*)))       ;raw spectroheliogram centre line
sphg=transpose(reform(cube(round(x0+offs),*,*)))  ;raw spectroheliogram continuum

if line_name eq 'HA' then begin
  sphe1=transpose(reform(cube(round(x0+2),*,*)))  ;raw spectroheliogram Ha-0.5A
  sphe2=transpose(reform(cube(round(x0-2),*,*)))  ;raw spectroheliogram Ha+0.5A
endif else begin
  sphe1=0
  sphe2=0
endelse

; calculation of rectified spectroheliogram

wdelete,2
DISKCOR3,sphg,spgcor,xdiamc,ydiamc              ; for continuum
DISKCOR3,sphe,sphcor,xdiam,ydiam,shft,xc,yc     ; for line centre

if line_name eq 'HA' then begin
  sph1cor=DSKC(sphe1,xdiam,ydiam,shft,xc,yc)    ; Ha wings rectification function,
  sph2cor=DSKC(sphe2,xdiam,ydiam,shft,xc,yc)    ;  using line centre parameters
endif else begin
  sph1cor=0
  sph2cor=0
endelse

; Dopplergram calculation

if line_name eq 'HA' then DIG_DOPPL,cube,x0,xdiam,ydiam,shft,xc,yc,doppl  $
  else begin
    doppl=0
    cube=0
  endelse

;------------------------------------------------------------------------
skok:

; horizontal strips and intensity asymmetry correction

if strip eq 'YES' then begin
    if line_name eq 'K3' then begin
       STRIPCOR,spgcor,corspg            ;K1 full correction
       STRIPCOR,sphcor,corsph,/skipstrip ;K3 only intensity asymmetry correction
    endif else begin
       STRIPCOR,spgcor,corspg            ;CN, HA full corrections
       STRIPCOR,sphcor,corsph
       STRIPCOR,sph1cor,cor1sph          ;HA wings
       STRIPCOR,sph2cor,cor2sph
    endelse
endif
if strip eq 'K3' then begin
    STRIPCOR,spgcor,corspg               ;full correction for all, also K3
    STRIPCOR,sphcor,corsph
    if line_name eq 'HA' then begin
       STRIPCOR,sph1cor,cor1sph          ;HA wings
       STRIPCOR,sph2cor,cor2sph
    endif
endif
if strip eq 'NO' then begin
    STRIPCOR,spgcor,corspg,/skipstrip    ;all only intensity asymmetry correction
    STRIPCOR,sphcor,corsph,/skipstrip    ;   corspg, corsph etc. are floats
    if line_name eq 'HA' then begin
       STRIPCOR,sph1cor,cor1sph,/skipstrip         ;HA wings
       STRIPCOR,sph2cor,cor2sph,/skipstrip
    endif

endif

; rotation of spectroheliograms to the position N-up, E-left.
;         drot - rotation angle (degrees counterclockwise)

spgr=rot(corspg,drot,cubic=(-0.5))
spgr=fix(spgr)>0   ; back to integer
sphr=rot(corsph,drot,cubic=(-0.5))
sphr=fix(sphr)>0   ; back to integer

if line_name eq 'HA' then begin
  sph1=rot(cor1sph,drot,cubic=(-0.5))
  sph1=fix(sph1)>0
  sph2=rot(cor2sph,drot,cubic=(-0.5))
  sph2=fix(sph2)>0
	dopp=rot(doppl, drot,cubic=(-0.5))	; float, in km/s
	dopp=fix(dopp*1000) 	; back to integer, !![ m/s ]!!
endif

; output ---------------------------------------------------

mlen=strlen(mereni)                    ;mereni contains the full path
measurement=strmid(mereni,mlen-9,8)    ;name of the measurement subdirectory
newdir=pathout+measurement

;Make a new output directory. If it already exists, destroy it and make a new one.

if file_test(newdir,/directory) eq 1 then SPAWN,'rmdir /s /q '+newdir
SPAWN,'mkdir '+newdir
outpath=newdir+'\'

wdelete,0
wset,3
xyouts,60,60,'Computing output images',/device,charsiz=2
wait,0.1

;FITS output

DIG_FITS4,sphr,outpath,measurement,line_name,rok,mesic,den,tfirst,tlast,drot,P0,ydiam,speed,progver
     ;centre line
DIG_FITS4,spgr,outpath,measurement,line_nam2,rok,mesic,den,tfirst,tlast,drot,P0,ydiamc,speed,progver
     ;continuum
if line_name eq 'HA' then begin
 DIG_FITS4,sph1,outpath,measurement,'HABLU',rok,mesic,den,tfirst,tlast,drot,P0,ydiamc,speed,progver
     ;Ha-0.5A
 DIG_FITS4,sph2,outpath,measurement,'HARED',rok,mesic,den,tfirst,tlast,drot,P0,ydiamc,speed,progver
     ;Ha+0.5A
 DOP_FITS,dopp,outpath,measurement,'DOP',rok,mesic,den,tfirst,tlast,drot,P0,ydiam,speed,progver
     ;dopplergram
endif

;JPEG output

DIG_JPEG,sphr,outpath,measurement,line_name,rok,mesic,den,tfirst,picli
     ;centre line
DIG_JPEG,spgr,outpath,measurement,line_nam2,rok,mesic,den,tfirst,picnt
     ;continuum
if line_name eq 'HA' then begin
 DIG_JPEG,sph1,outpath,measurement,'HABLU',rok,mesic,den,tfirst,pichb
 DIG_JPEG,sph2,outpath,measurement,'HARED',rok,mesic,den,tfirst,pichr
 DOP_JPEG,dopp,outpath,measurement,'DOP',rok,mesic,den,tfirst,picdop
     ;dopplergram
endif

;display on main screen

wset,3
aspect=1000./1200.
sx=!d.x_size/5    ;x-size of pictures to be displayed in window 3
sy=fix(sx*aspect)

if line_name eq 'K3' then begin

  for i=0,4 do tvscl,congrid(picli(*,*,i),sx,sy),i
  for i=0,4 do tvscl,congrid(picnt(*,*,i),sx,sy),i+5
  tvscl,congrid(picli(*,*,5),sx,sy),14
endif

if line_name eq 'HA' then begin
 picha=bytarr(1200,1000,10)
 picha(*,*,0:1)=pichb(*,*,1:2)   ;2x HABL
 picha(*,*,2:4)=picli(*,*,1:3)   ;3x HA
 picha(*,*,5:6)=pichr(*,*,1:2)   ;2x HARE
 picha(*,*,7:8)=picnt(*,*,1:2)   ;2x CN
 picha(*,*,9)=picli(*,*,5)       ; Z
 for i=0,9 do tvscl,congrid(picha(*,*,i),sx,sy),i
 tvscl,congrid(picdop,sx,sy),14
 picli=picha(*,*,0:4)            ;arrays for zooming
 picnt=picha(*,*,5:9)
endif
xyouts,60,60,'Computing output images',/device,charsiz=2,color=0

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

if runmo eq 1 or runmo eq 3 then begin

   xlim=lonarr(6)
   ylim=lonarr(3)
   for i=0,5 do xlim(i)=sx*i
   for i=0,2 do ylim(i)=!d.y_size-(i+1)*sy
   xyouts,60,160,'Click images to zoom or empty area to continue',/device,charsiz=3,charthick=3

   curs:

   cursor,x,y,/device,/down
   if y lt ylim(2) then goto,fin
   if y lt ylim(1) and x lt xlim(4) then goto,fin
   if y ge ylim(0) then begin
     i=(-1)
     repeat i=i+1 until (x ge xlim(i)) and (x lt xlim(i+1))
     DIG_TVW,picli(*,*,i)
     goto,curs
   endif
   if (y lt ylim(0)) and (y ge ylim(1)) then begin
     i=(-1)
     repeat i=i+1 until (x ge xlim(i)) and (x lt xlim(i+1))
     DIG_TVW,picnt(*,*,i)
   endif else begin
	if line_name eq 'HA' then DIG_TVW,picdop else DIG_TVW,picli(*,*,5)
   endelse
   goto,curs

endif else wait,5

fin:
END