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