! @(#)pisco.prg	17.1.1.1 (ES0-DMD) 01/25/02 17:55:20
!*********************************************************************
!*
!*
!*                            PISCO
!*                            -----
!*
!*  Purpose:   reduces polarimetric data obtained with PISCO
!*
!*   author:   Otmar Stahl            March 1986
!*              updated to portable MIDAS by
!*             Martin Schloetelburg    June 1990
!*
!*      use:   @@ PISCO catalog result sky calibration mode
!*             
!*             catalog     = catalog containing data files to be reduced
!*             result      = output table
!*             sky         = sky measurement file
!*             calibration = sky measurement file
!*             mode     = 1 : x and y channel are reduced seperately
!*                        2 :     ''          ''          together
!*                        3 :     ''          ''          both (1+2)
!*
!*  exampel:   @@ PISCO CATALOG TABLE SKY CALIB 3
!*
!*             This would reduce all the files listed in catalog CATALOG
!*             taking the file SKY as sky measurement and CALIB for calibration
!*             with x and y channel both separately and together. The result 
!*             is stored in the table file TABLE.
!*
!* defaults:   catalog = ? , result = ? , sky = ? , calibration = ? , mode = 3
!*
!*     note:   The error is estimated from the noise level in the 
!*             Fourier-transform. The angle is still arbitrary and has 
!*             to be transformed to a standard system by observations of 
!*             standard stars.
!*
!*             All commands are explained in the MIDAS manual.
!*             Scratch files are called XX* and are deleted afterwards.
!*
!*
!*********************************************************************
!
DEFINE/PARA P1 ? ? "catalog of input data files"
DEFINE/PARA P2 ? ? "output table holding results"
DEFINE/PARA P3 ? ? "input sky file"
DEFINE/PARA P4 ? ? "input calibration file"
DEFINE/PARA P5 3 ? "data reduction mode"
!----------------------------------
!
! get number of data files and prevent 
! further entries in image catalog
!
!----------------------------------
SHOW/ICAT P1 NODISPLAY
CLEAR/ICAT
!----------------------------------
!
! create output table
!
!----------------------------------
COMP/KEYW  OUT_A(1:20) = M$TIME()
WRITE/OUT {OUT_A(13:20)} creating output table {P2}
WRITE/OUT
!
CREA/TABLE  {P2} 17 {OUTPUTI(1)} NULL
CREA/COLUMN {P2} :IDENT "A" A16 C*16
CREA/COLUMN {P2} :UT "A" E12.5
CREA/COLUMN {P2} :X_intensity "A" E12.5
CREA/COLUMN {P2} :Y_intensity "A" E12.5
CREA/COLUMN {P2} :polxy "A" E12.5
CREA/COLUMN {P2} :errorxy "A" E12.5
CREA/COLUMN {P2} :anglexy "A" E12.5
CREA/COLUMN {P2} :polx "A" E12.5
CREA/COLUMN {P2} :errorx "A" E12.5
CREA/COLUMN {P2} :anglex "A" E12.5
CREA/COLUMN {P2} :poly "A" E12.5
CREA/COLUMN {P2} :errory "A" E12.5
CREA/COLUMN {P2} :angley "A" E12.5
!-------------------------------------------------
!
! extract useful part of sky
!
!-------------------------------------------------
EXTRA/LIN XXS1 = {P3}[1,1:32,1]
EXTRA/LIN XXS2 = {P3}[1,2:32,2]
!-------------------------------------------------
!
! extract normalized calibration
!
!-------------------------------------------------
EXTRA/LIN XXSC = {P4}[1,8:32,8]
!-------------------------------------------------
!
! divide sky by integration time
!
!-------------------------------------------------
COMP/IMAG XXS1 = XXS1/{{P3},O_TIME(7)}
COMP/IMAG XXS2 = XXS2/{{P3},O_TIME(7)}
!-------------------------------------------------
!
! do-loop for object files
!
!-------------------------------------------------
DEFINE/LOCAL CATAL/I/1/1 0
CAT_LOOP:
!
! get name of data files from catalog
!
STORE/FRAME IN_A {P1}.cat 1 FINITO
!
COMP/KEYW  OUT_A(1:20) = M$TIME()
WRITE/OUT {OUT_A(13:20)} processing file {IN_A}
!
! store file identifier in table column :IDENT
!
{P2},:IDENT,@{CATAL} = "{{IN_A},IDENT(1:16)}"
!
! store UT in table
!
{P2},:UT,@{CATAL} = {{IN_A},O_TIME(5)}
!
! extract useful part
!
EXTRA/LIN XXI1 = {IN_A}[1,1:32,1]
EXTRA/LIN XXI2 = {IN_A}[1,2:32,2]
!
! divide by integration time, subtract sky and calibrate
!
COMP/IMAG XXI1 = (XXI1/{XXI1,O_TIME(7)}-XXS1)/XXSC
COMP/IMAG XXI2 = (XXI2/{XXI1,O_TIME(7)}-XXS2)/XXSC
!
! normalize and store mean value
!
STAT/IMAG XXI1 ? ? ? RN
COMP/IMAG XXI1 = XXI1/{OUTPUTR(3)}-1.
{P2},:X_intensity,@{CATAL} = {OUTPUTR(3)}
!
STAT/IMAG XXI2 ? ? ? RN
COMP/IMAG XXI2 = XXI2/{OUTPUTR(3)}-1.
{P2},:Y_intensity,@{CATAL} = {OUTPUTR(3)}
!
! branch to modes
!
IF P5 .EQ. 2 THEN
  GOTO XANDY
ENDIF
!
! reduce X channel data
!
FFT/FPOW  XXI1
COMP/KEYW OUTPUTR(1) = 205.234*{POWER[@13]}
{P2},:polx,@{CATAL} = {OUTPUTR(1)}
!
! compute error
!
COMP/IMAG XXIE = POWER*POWER
STAT/IMAG XXIE ? ? ? RN
COMP/KEYW OUTPUTR(1) = M$SQRT({OUTPUTR(3)})*232.2
{P2},:errorx,@{CATAL} = {OUTPUTR(1)}
!
! compute angle
!
IF FFTR[@13] .GT. 0. THEN
COMP/KEYW OUTPUTR(1) = 45.-0.5*M$ATAN({FFTI[@13]}/{FFTR[@13]})
{P2},:anglex,@{CATAL} = {OUTPUTR(1)}
ELSE
COMP/KEYW OUTPUTR(1) = 135.-0.5*M$ATAN({FFTI[@13]}/{FFTR[@13]})
{P2},:anglex,@{CATAL} = {OUTPUTR(1)}
ENDIF
!
! reduce Y channel data
!
FFT/FPOW  XXI2
COMP/KEYW OUTPUTR(1) = 205.234*{POWER[@13]}
{P2},:poly,@{CATAL} = {OUTPUTR(1)}
!
! compute error
!
COMP/IMAG XXIE = POWER*POWER
STAT/IMAG XXIE ? ? ? RN
COMP/KEYW OUTPUTR(1) = M$SQRT({OUTPUTR(3)})*232.2
{P2},:errory,@{CATAL} = {OUTPUTR(1)}
!
! compute angle
!
IF FFTR[@13] .LT. 0. THEN
COMP/KEYW OUTPUTR(1) = 45.-0.5*M$ATAN({FFTI[@13]}/{FFTR[@13]})
{P2},:angley,@{CATAL} = {OUTPUTR(1)}
ELSE
COMP/KEYW OUTPUTR(1) = 135.-0.5*M$ATAN({FFTI[@13]}/{FFTR[@13]})
{P2},:angley,@{CATAL} = {OUTPUTR(1)}
ENDIF
!
! now use differential method
!
IF P5 .EQ. 1 THEN
GOTO NEXT
ENDIF
XANDY:
!
COMP/IMAG XXI3 = (XXI1-XXI2)/(XXI1+XXI2+2)
!
! Fourier-transform of two-channel data
!
FFT/FPOW  XXI3
COMP/KEYW OUTPUTR(1) = 205.234*{POWER[@13]}
{P2},:polxy,@{CATAL} = {OUTPUTR(1)}
!
! compute error
!
COMP/IMAG XXIE = POWER*POWER
STAT/IMAG XXIE ? ? ? RN
COMP/KEYW OUTPUTR(1) = M$SQRT({OUTPUTR(3)})*232.2
{P2},:errorxy,@{CATAL} = {OUTPUTR(1)}
!
! compute angle
!
IF FFTR[@13] .GT. 0. THEN
COMP/KEYW OUTPUTR(1) = 45.-0.5*M$ATAN({FFTI[@13]}/{FFTR[@13]})
{P2},:anglexy,@{CATAL} = {OUTPUTR(1)}
ELSE
COMP/KEYW OUTPUTR(1) = 135.-0.5*M$ATAN({FFTI[@13]}/{FFTR[@13]})
{P2},:anglexy,@{CATAL} = {OUTPUTR(1)}
ENDIF
!
NEXT:
!
DELETE XXI%  NO
DELETE FFTI  NO
DELETE FFTR  NO
DELETE POWER NO
!
GOTO CAT_LOOP
FINITO:
DELETE XXS% NO
RETURN