;---------------------------------------------------------------------------
; Document name: polarization.pro
; Created by:    Oliver Trachsel, January 2002
;
; Last Modified: 2002/01/30
;---------------------------------------------------------------------------
;
;+
; NAME:
;       Automatic polarization determination with background subtraction
;
; PURPOSE:
;       Determine the L/R-Polarization of a data point using the
;       following steps:
;		1) Let user choose a point and two reference bands
;                  for background subtraction
;		2) Search and open the p- and i-file(s) containing the
;                  specified time/frequency data
;		3) Calculate R and L intensities for data point and
;                  background reference bands, subtract background,
;                  calculate new polarization
;
; CALLING SEQUENCE
;	DriftRate, zReal, xAxis, yAxis, WINNR=WINNR
;
; INPUTS:
;       xAxis: the x-axis of the displayed image zone
;       yAxis: the corresponding y-axis 
;
; KEYWORDS:
;       WINNR: If set, the specified window contains the image to be
;              treaten, otherwise the active window is taken.
;
; SIDE EFFECTS:
;
; RESTRICTIONS:
;
; MODIFICATION HISTORY:
;       02/01/30 Initial version
;                Oliver Trachsel <toliver@astro.phys.ethz.ch>
;		 Institute of Astronomy, ETH Zurich
;-

; utility procedure to swap contents of x and y
PRO SwapVars, x, y
    tmp = x
    x = y
    y = tmp
END

PRO Polarization, xAxis, yAxis, WINNR=winNr

    COMMON ImageReal, img,xglobl,yglobl, fName ; Data

    IF N_ELEMENTS(WINNR) NE 1 THEN winNr=!D.WINDOW

; determine the date of the current plot, first look at
; the beginning of the filename, otherwise ask user

    year = STRMID(fName,0,4)
    month = STRMID(fName,4,2)
    day = STRMID(fname,6,2)
    year_num = LONG(year)
    month_num = LONG(month)
    day_num = LONG(day)

    WHILE (year_num LT 1900) OR (month_num LT 1) OR (month_num GT 12) OR $
      (day_num LT 1) OR (day_num GT 31) DO BEGIN

        PRINT
        datestring=''
        READ, datestring, $
          PROMPT="Please enter the date of the current data (YYYY/MM/DD): "
        PRINT
        year = STRMID(datestring,0,4)
        month = STRMID(datestring,5,2)
        day = STRMID(datestring,8,2)
        year_num = LONG(year)
        month_num = LONG(month)
        day_num = LONG(day)
        IF STRLEN(datestring) NE 10 THEN year_num = -1
    ENDWHILE
    date = year+'/'+month+'/'+day

    GetImageSizes, xPos, yPos, width, height
    GetOffsets, d1, imageOffset
    GetBreaks, dxBreaks, dyBreaks    
    nx = N_Elements(xAxis)
    ny = N_Elements(yAxis)

; let user choose data point and two background intervals

ChooseCoords:    
    XMessage, ['- Left mouse button to choose data point', $
               '- Right button to exit'], $
      TITLE='Choose point', XOFF=0, YOFF=0, WBASE=wbase
    WSHOW, winNr
    WSET, winNr
    Cursor, xPt, yPt, /DEVICE, /UP
    XKill, wbase
    IF !MOUSE.BUTTON EQ 4 THEN RETURN
    
    XMessage, ['Choose start point for first', $
               'background interval'], $
      TITLE='Choose point', XOFF=0, YOFF=0, WBASE=wbase
    WSHOW, winNr
    WSET, winNr
    Cursor, x0Bg1, y, /DEVICE, /UP
    XKill, wbase
    
    XMessage, ['Choose end point for first', $
               'background interval'], $
      TITLE='Choose point', XOFF=0, YOFF=0, WBASE=wbase
    WSHOW, winNr
    WSET, winNr
    Cursor, x1Bg1, y, /DEVICE, /UP
    XKill, wbase
    
    XMessage, ['Choose start point for second', $
               'background interval'], $
      TITLE='Choose point', XOFF=0, YOFF=0, WBASE=wbase
    WSHOW, winNr
    WSET, winNr
    Cursor, x0Bg2, y, /DEVICE, /UP
    XKill, wbase
    
    XMessage, ['Choose end point for second', $
               'background interval'], $
      TITLE='Choose point', XOFF=0, YOFF=0, WBASE=wbase
    WSHOW, winNr
    WSET, winNr
    Cursor, x1Bg2, y, /DEVICE, /UP
    XKill, wbase
        
; convert the clicked coordinates to data coordinates
    
    IF x0Bg1 GT x1Bg1 THEN SwapVars, x0Bg1, x1Bg1

    IF x0Bg2 GT x1Bg2 THEN SwapVars, x0Bg2, x1Bg2
    IF (x1Bg1 GE x0Bg2) AND (x1Bg2 GE x0Bg1) THEN BEGIN
        PRINT, 'Overlapping backround intervals. Restarting.'
        GOTO,ChooseCoords
    ENDIF
    
    xPt = DevToPix(xPt, xPos, width, dxBreaks, imageOffset, nx)
    yPt = DevToPix(yPt, yPos, height, dyBreaks, imageOffset, ny)
    x0Bg1 = DevToPix(x0Bg1, xPos, width, dxBreaks, imageOffset, nx)
    x1Bg1 = DevToPix(x1Bg1, xPos, width, dxBreaks, imageOffset, nx)
    x0Bg2 = DevToPix(x0Bg2, xPos, width, dxBreaks, imageOffset, nx)
    x1Bg2 = DevToPix(x1Bg2, xPos, width, dxBreaks, imageOffset, nx)
; convert to data coordinates
    xPt_val = xAxis[xPt]
    yPt_val = yAxis[yPt]
    x0Bg1_val = xAxis[x0Bg1]
    x1Bg1_val = xAxis[x1Bg1]
    x0Bg2_val = xAxis[x0Bg2]
    x1Bg2_val = xAxis[x1Bg2]

; try to open files for the whole active zone and store values in a
; local variable
    
    IF N_ELEMENTS(zoneRead) LE 0 THEN BEGIN
        PRINT,'Reading intensity and polarization data for active area...'
        view_time_intv = date + ' ' + [HMSConvert(xAxis[0]), $
                                       HMSConvert(xAxis[nx-1])]
        view_freq_intv = [yAxis[ny-1], yAxis[0]]
        view_pol = rapp_get_spectrogram(view_time_intv, view_freq_intv, /POL, $
                                  XAXIS = pol_xAxis, YAXIS = pol_yAxis)
        IF NOT (DATATYPE(view_pol) EQ 'STR') THEN $
          view_int = rapp_get_spectrogram(view_time_intv, view_freq_intv, $
                                          XAXIS = int_xAxis, YAXIS = int_yAxis)
        zoneRead = NOT ((DATATYPE(view_pol) EQ 'STR') OR $
                        (DATATYPE(view_int) EQ 'STR'))
        IF zoneRead THEN PRINT, '  ... success!' ELSE PRINT, '  ... failed!'
        PRINT
    END

; extract the necessary intervals
    IF zoneRead THEN BEGIN
        pt_i = view_int[xPt, yPt]
        bg1_i = view_int[x0Bg1:x1Bg1, yPt]
        bg2_i = view_int[x0Bg2:x1Bg2, yPt]
        pt_p = view_pol[xPt, yPt]
        bg1_p = view_pol[x0Bg1:x1Bg1, yPt]
        bg2_p = view_pol[x0Bg2:x1Bg2, yPt]
        
; if there wasn't data available for the whole displayed zone
; then try to read only the specified point and intervals
    ENDIF ELSE BEGIN    ; i.e. if NOT zoneRead
        ; create intervals 
        time_intv_Pt = date+' '+[HMSConvert(xPt_val),HMSConvert(xPt_val)]
        time_intv_Bg1 = date+' '+[HMSConvert(x0Bg1_val),HMSConvert(x1Bg1_val)]
        time_intv_Bg2 = date+' '+[HMSConvert(x0Bg2_val),HMSConvert(x1Bg2_val)]
        freq_intv = [yPt_val, yPt_val]

        ; try to get all necessary data, if not ok then restart
        bg1_p = rapp_get_spectrogram(time_intv_Bg1, freq_intv, /POL, $
                                     XAXIS=bg1_xaxis, YAXIS=bg1_yaxis)
        bg2_p = rapp_get_spectrogram(time_intv_Bg2, freq_intv, /POL, $
                                     XAXIS=bg2_xaxis, YAXIS=bg2_yaxis)
        pt_p = rapp_get_spectrogram(time_intv_Pt, freq_intv, /POL, $
                                    XAXIS=pt_xaxis, YAXIS=pt_yaxis)
        bg1_i = rapp_get_spectrogram(time_intv_Bg1, freq_intv, $
                                     XAXIS=bg1_xaxis, YAXIS=bg1_yaxis)
        bg2_i = rapp_get_spectrogram(time_intv_Bg2, freq_intv, $
                                     XAXIS=bg2_xaxis, YAXIS=bg2_yaxis)
        pt_i = rapp_get_spectrogram(time_intv_Pt, freq_intv, $
                                    XAXIS=pt_xaxis, YAXIS=pt_yaxis)
        
        IF (DATATYPE(pt_p) EQ 'STR') OR (DATATYPE(bg1_p) EQ 'STR')  $
          OR (DATATYPE(bg2_p) EQ 'STR') OR (DATATYPE(pt_i) EQ 'STR') $
          OR (DATATYPE(bg1_i) EQ 'STR') OR (DATATYPE(bg2_i) EQ 'STR') $
          THEN BEGIN
            PRINT, 'Error: Not all necessary data available. Restarting'
            GOTO,ChooseCoords
        ENDIF
    ENDELSE

; do the calculations (extraction of L and R,
; background subtraction)

    i = [bg1_i, bg2_i, pt_i] ; rapp_get_spectrogram returns linearized i
    p = [bg1_p, bg2_p, pt_p]
    tmp = (100.0-p)/(100.0+p)
    denominator = 1.0 + tmp
    R = i/denominator
    L = R*tmp

    bg1_max = N_ELEMENTS(bg1_i)-1
    bg2_min = bg1_max+1
    bg2_max = N_ELEMENTS(i)-2
    
    R_bg1 = (N_ELEMENTS(bg1_i) GT 1) ? MEAN(R[0:bg1_max], /DOUBLE) $
      : R[0]
    L_bg1 = (N_ELEMENTS(bg1_i) GT 1) ? MEAN(L[0:bg1_max], /DOUBLE) $
      : L[0]
    R_bg2 = (N_ELEMENTS(bg2_i) GT 1) ? MEAN(R[bg2_min:bg2_max], /DOUBLE) $
      : R[bg2_min]
    L_bg2 = (N_ELEMENTS(bg2_i) GT 1) ? MEAN(L[bg2_min:bg2_max], /DOUBLE) $
      : L[bg2_min]

    ; interpolate (or extrapolate) background to x-coord of chosen data point
    xmid_bg1 = x0Bg1+(x1Bg1-x0Bg1)/2.0
    xmid_bg2 = x0Bg2+(x1Bg2-x0Bg2)/2.0
    R_slope = (R_bg2 - R_bg1)/(xmid_bg2-xmid_bg1)
    L_slope = (L_bg2 - L_bg1)/(xmid_bg2-xmid_bg1)
    Rbg = R_slope*(xPt-xmid_bg1)+R_bg1
    Lbg = L_slope*(xPt-xmid_bg1)+L_bg1
    
;    old version without interpolation
;    Rbg = MEAN(R[0:N_ELEMENTS(R)-2], /DOUBLE)
;    Lbg = MEAN(L[0:N_ELEMENTS(L)-2], /DOUBLE)

    Rsig = R[N_ELEMENTS(R)-1] - Rbg
    Lsig = L[N_ELEMENTS(L)-1] - Lbg
    p_sig = 100.0*((Rsig-Lsig)/(Rsig+Lsig))

; output the results
    PRINT, '--------------------------------------------------------------------'
    PRINT, '  POLARIZATION DETERMINATION'
    PRINT, '--------------------------------------------------------------------'
    PRINT
    PRINT, '  Date: ',date
    PRINT
    PRINT, '  Chosen coordinates:'
    PRINT, '    Data point:   ',HMSConvert(xPt_val),' / ',$
      STRTRIM(STRING(yPt_val, FORMAT='(D10.1)'),1),' MHz'
    PRINT, '    Background 1: ['+HMSConvert(x0Bg1_val)+', '+HMSConvert(x1Bg1_val)+']'
    PRINT, '    Background 2: ['+HMSConvert(x0Bg2_val)+', '+HMSConvert(x1Bg2_val)+']'
    PRINT
    PRINT, '  Results:'
    PRINT, '    Interpolated background: L = ',STRN(Lbg),' sfu, R = ',STRN(Rbg),' sfu'
    PRINT, '    Pixel value after subtr: L = ',STRN(Lsig),' sfu, R = ', $
      STRN(Rsig),' sfu'
    PRINT
    PRINT, '    Polarization: ',STRN(p_sig),' %'
    PRINT
    PRINT, '--------------------------------------------------------------------'
    PRINT
    GOTO,ChooseCoords
END