! 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 ----------------"
*