Listing of the stokes.class procedure

!    stokes.class
!
! written by : P. Colom, with help from Eric Gérard
!              Feb. 2008 
! aims: stokes parameters, linear polarisation parameter, PA, (U/Q)
!
! you need 2 scans: first with 0 and 90 deg   (Rmer = -45)
!		    second    +45    -45      (Rmer =   0)
!
! to fit a baseline, you need to mask the line and the spectrum start/end 
! channels
define character outfile*20, infile*20 
define integer first_scan last_scan
define real Vel_start Vel_end Vmin Vmax
!
GREG1\SET /DEFAULT
GREG1\PENCIL /DEFAULT
say "input file (MAX 20 char) ? "
let infile = 
say "output file (MAX 20 char, new name) ? "
let outfile = 
!
file in 'infile'
file out 'outfile' new
!
say "first scan (0°,90°) :"
let first_scan = 
say "last scan (+45°,-45°) :"
let last_scan =
!
set scan first_scan first_scan       ! angles : first, second = 0, 90°
find /all
list
set format brief    ! brief header for plot
set weigh equal     ! equal weight for sum or accu
get f
! define dimension of angle psi
define real psi[channels]
plot
get n
plot
accu                ! Stokes I
plot
say " eliminate spectrum start/end channels : Vel_start and Vel_end"
let Vel_start = 
let Vel_end = 
say " baseline fit : you need to mask the line between Vmin & Vmax"
let Vmin = 
let Vmax = 
!
set mode x Vel_start Vel_end       ! eliminates channels at both ends
set window Vmin Vmax               ! masks the line
base 2 /plot                   ! second order polynomial fit
! pause " Stokes I -- type cont to go ahead"
!
memorize I1
get f
get n
multiply -1
accu 		    ! Stokes Q
plot
base 2 /plot
memorize SQ
! pause " Stokes Q -- type cont to go ahead"
get n 		    ! LCP
get n		    ! RCP
multiply -1
accu                ! Stokes V
plot 
base 2 /plot
memorize V1
! pause " Stokes V -- type cont to go ahead"
!  and now: rotation of feed horn by 45 deg
set scan last_scan last_scan       ! angles : first, second = = +45, -45
find /all
list 
get f
plot 
get n
plot
accu                ! Stokes I
plot
base 2 /plot
memorize I2
! pause " Stokes I -- type cont to go ahead"
get f
get n
multiply -1
accu                ! Stokes U
plot 
base 2 /plot
memorize SU
! pause " Stokes U -- type cont to go ahead"
get n               ! LCP
get n               ! RCP
multiply -1
accu                ! Stokes V
plot 
base 2 /plot
memorize V2
! pause " Stokes V -- type cont to go ahead"
!
! compute averages of I and V, and ratio U/Q
retrieve I1
retrieve I2
accu
mult 0.5
memorize SI         ! Stokes I
plot 
write
pause " average of Stokes I -- type cont to go ahead"
retrieve SQ
write
plot
pause " Stokes Q -- type cont to go ahead"
!  degree of linear polarization
memorize Pl         ! Q in Pl
retrieve Pl
let RY = RY*RY      ! Q^2
memorize Pl
retrieve SU
let RY = RY*RY      ! U^2
retrieve Pl
accu                ! Q^2 + U^2
let RY = sqrt(RY)
retrieve SI
swap                !  Exchange the contents of the R and T buffers
divide 0.5          !  sqrt(Q^2 + U^2) / I
memorize Pl
!
retrieve SQ
retrieve SU
write
plot 
pause " Stokes U -- type cont to go ahead"
divide 0.5           ! divide U by Q  (0.5 is a threshold, avoid /0)
plot
pause " U/Q  -- type cont to go ahead"
!
psi = 0.5*atan(ry)        ! position angle
psi = psi*180./pi
ry = psi
set mode y -90 90
plot
write
pause " psi = 1/2 arctan(U/Q) -- type cont to go ahead"
!
retrieve Pl
set mode y 0 1
plot
write
pause " degree of linear polarization -- type cont to go ahead"
!
retrieve V1
retrieve V2
accu
mult 0.5
set mode y total
plot
write
pause " average of Stokes V -- type cont to go ahead "
set scan * *            ! release constraint on scan number
say " ---------------------------------------------------"
say " conclusion: we have written in " 'outfile'
say " 1) "
say " 2)  Q "
say " 3)  U "
say " 4)  psi = 1/2 arctan(U/Q) "
say " 5)  degree of linear polarization"
say " 6) "
say " -------------- end of stokes.class ----------------"