poltools: module of polarimetric tools¶
PyHdust poltools module: polarimetry tools
History: grafpol working for *_WP1110….log files! grafpol working for log/out files with more than a single star
coauthor:  Daniel Bednarski 

license:  GNU GPL v3.0 https://github.com/danmoser/pyhdust/blob/master/LICENSE 

pyhdust.poltools.
chkStdLog
(f, calc, path=None, delta=3.5, verbose=True)[source]¶ Verify if there are standards for filter f and calcite calc inside path/std.dat. Return True if successful, unless the standard has been reduced and marked with E flag.
delta is the allowed variation for the angles between the two beams for one same calcite.

pyhdust.poltools.
chooseout
(objdir, obj, f, nstar=1, sigtol=<function <lambda>>)[source]¶ Olha na noite, qual(is) *.OUT(s) de um filtro que tem o menor erro.
Retorna um mais valores para toda a sequencia (i.e., pasta). Criterios definidos no anexo de polarimetria.
minerror == True: recebe o out de menor erro em qualquer caso.
sigtol: condicao de tolerancia para pegar o agrupamento de menor erro. Se o sigma do agrupamento com todas N posicoes for menor que a funcao sigtol sobre o sigma do agrupamento de menor erro, entao usa o agrupamento com as N posicoes; do contrario usa o de menor erro. Exemplo com 16 posicoes: erro do grupo de 16 posicoes == 0.000230;
menor erro == 0.000200. Se sigtol(sig) = 1.4*sig, como sigtol(0.000200) == 0.000240 > 0.000230, usa o agrupamento de 16 posicoes.O numero de posicoes eh baseado no numero de arquivos *.fits daquele filtro.
 ‘nstar’ == star number inside ‘out’ file (usefull when there are
 more than a single star inside .out)

pyhdust.poltools.
corObjStd
(night, f, calc, path=None, delta=3.5, verbose=True)[source]¶ Find the correction factor delta theta for filter ‘f’ inside night ‘night’, for calcite ‘calc’ (the last variable must be the angle between the ord. and extraord. beams). Returns the values for matching standard stars inside ‘night’/std.dat file, except when marked with an ‘E’ flag.
 NEW: In case of missing data in filter ‘f’, this routine
 tries to compute delta theta by using the data from another filter, making use of the relations among the values of delta theta for each filter). But only the estimate presenting the smaller error will be taken.
 verbose: auxiliary variable used as False inside
 polt.genTarget routine. Keep it as True otherwise.
 Formulas:
dth = fact*th^std_measured  th^std_published
 fact = 1 for observations before 2015 March 1st
+1 otherwise
 delta: tolerance, in degree, of angle of two beams for
 one same calcite (default: +/3.5 degree)
 Output, in order: stdnames, mdth, smdth, flag, tags
 stdnames: string with standard star names separated
 by commas (,)
 mdth: mean delta theta (corrected by the output factor
 +1 of routine polt.thtFactor())
 smdth: the error of the mean delta theta
 flag: flag concerning to the standards.
 tags: string with the tags concerning to the standards,
 separated by commas (,), or ‘—’ if none.

pyhdust.poltools.
countStars
(objdir, f)[source]¶ Count how many stars there are inside outfiles in ‘objdir’ and filter f. Return 0 if there are no outfiles

pyhdust.poltools.
fitMCMCline
(x, y, sx, sy, star='', margin=False, plot_adj=True, fig=None, ax=None, n_burnin=350, n_mcmc=600, n_walkers=120, thet_ran=[0.0, 180.0], b_ran=[1.0, 1.0], Pb_ran=[0.0, 1.0], Yb_ran=[1.0, 1.0], Vb_ran=[0.0, 1.0], extens='pdf')[source]¶ Fit a line using Markov Chain Monte Carlo for data with both x and y errors and with bad points from emcee code.
The model is tha sum of two models: a) a line model for the good points, y = ax + b, with data displaced ortogonally by a gaussian diplacentment (covariance is suposed null!); b) a gaussian distribution orthogonal to the line for the bad points with amplitude, mean and variance equal to Pb, Yb and Vb (see Hoog, Bovy and Lang, ``Data analysis recipes: Fitting a model to data’‘)
The MCMC does a ``mixture’’ among both model for each point. So, it is not needed know about the bad and good points necessairly! The five parameters to be found are theta=arctan(a), b, Pb, Yb and Vb.
INPUT:
 x/y/sx/sy: array/list with the data
 star: star name to be printed in the graph and
 its filename. If it’s a void str ‘’, this routine give a random number to prevent overwriting data.
 margin: marginalize the corner graphs over Pb,
 Yb, Vb (generating the graph only for theta and b)?
 plot_adj: show a plot of data+fit?
 fig: Figure to append the plots for data+fit.
 If None, a new figure is generated.
ax: Axes, as like fig above.
 n_burnin: number of iterations for burningin
 n_mcmc: number of iterations to run emcee
 n_walkers: number of walkers to map the posterior
 probabilities.
 thet_ran: [thet_min, thet_max]
 b_ran: [b_min, b_max]
Pb_ran: [Pb_min, Pb_max], with Pb_min >= 0 Yb_ran: [Yb_min, Yb_max] Vb_ran: [Vb_min, Vb_max], with Vb_min >= 0 extens: extension for the graph file
OUTPUT:
 b*cos(theta): Idem, for b*cos(theta).
 Pb: Idem, for Pb. Yb: Idem, for Yb. Vb: Idem, for Vb.
 fig1: Figure pointer to the corner graph ([]
 if show==False).
 fig2: Figure pointer to the data+fit graph ([]
 if plot_adj==False).
FORMULAS:
Supposing i as each data point, the log of Likelihood function (L) takes form:
log(L) = sum(log(p_good_i+p_bad_i))with
 p_good_i = (1Pb)/sqrt(2*pi*var_i)*
 exp(0.5*(disp_i**2/var_i)),
 p_bad_i = Pb/sqrt(2*pi*(Vb+var)*
 exp(0.5*(disp_iYb)**2)/(Vb+var_i))
where disp_i is the total projection of the (x_i,y_i) values over the line and var_i, the projected variance:
 disp_i = v*Z_i b cos(thet)
 var_i = v*S_i*v
with
 v = (sin(theta), cos(theta)) (versor orthogonal
 to the line)
Z_i = (x_i, y_i) (data point) S_i =  sx_i^2 sxy_i^2 = sx_i^2 0  = (covariance
syx_i^2 sy_i^2  0 sy_i^2 matrix)

pyhdust.poltools.
fixISP
(logfile, ispol, path2=None)[source]¶ Correct the interstellar polarization in file logfile
logfile: outfile from genTarget. ispol: the Serkowski parameters [P_max, lambda_max, theta_IS]
to correct IS polarization (P_max ein % and lambda_max in Angstrom). If ispol==None, don’t make the correction of ISP. path2: path where to save the data file. If None, it is
 supposed the same of logfile.

pyhdust.poltools.
genAllLog
(path=None, sigtol=<function <lambda>>, autochoose=False, delta=3.5)[source]¶ Generate the std.dat/obj.dat for one reduced night
path: path of night delta: tolerance for the angle between the two beams of calcite.
If abs(angle1  angle2) < delta, both observations 1 and 2 are assigned as the same calcite. sigtol: tolerance to use the outfiles with all WP instead the
 out with best error. Must be a ‘function’ that receives a pol sigma value (in decimal system and NOT in per cent, i.e., value from 0. to 1., where 1 is a 100% polarized source) and return the maximum sigma for which to ignore the best out. Its format must be a ‘python’s lambda function’! The default values is sigtol=lambda sigm: 1.4*sigm, while the old value was sigtol=lambda sigm: 1.1*sigm + 0.00005. If you want to take just the groups with all WP, and none other, you can specify sigtol=lambda sigm: 1000.*sigm, for example.
 autochoose: choose best outfiles automatically, without
 interaction?

pyhdust.poltools.
genLog
(path, subdirs, tgts, fileout, sigtol=<function <lambda>>, autochoose=False, delta=3.5)[source]¶ Generate the .dat file with data of objects ‘tgts[:]’ inside ‘path’/’subdirs[:]’ directories
Save the results in ‘path’/’fileout’ Usable to generate target and standard lists (obj.dat and std.dat)
 delta: tolerance for the angle between the two beams of calcite.
 If abs(angle1  angle2) < delta, both observations 1 and 2 are assigned as the same calcite.
 sigtol: tolerance to use the outfiles with all WP instead the
 out with best error. Must be a ‘function’ that receives a pol sigma value (in decimal system and NOT in per cent, i.e., value from 0. to 1., where 1 is a 100% polarized source) and return the maximum sigma for which to ignore the best out. Its format must be a ‘python’s lambda function’! The default values is sigtol=lambda sigm: 1.4*sigm, while the old value was sigtol=lambda sigm: 1.1*sigm + 0.00005. If you want to take just the groups with all WP, and none other, you can specify sigtol=lambda sigm: 1000.*sigm, for example.
 autochoose: choose best outfiles automatically, without
 interaction?

pyhdust.poltools.
genTarget
(target, path=None, path2=None, ispol=None, skipdth=False, delta=3.5, epssig=2.0)[source]¶ Gen. target
Generate a table with all observations found for ‘target’, unless those which have E flag.
The error in delta theta (factor from observation of standard star) IS NOT propagated to the theta of Be star.
 skipdth: skip the observations without estimates for delta
 theta (no standard star in none filter and no standard star in previous and next nights)? Case True, the exception is only when the data is ‘unpolarized’ (defined by ‘epssig’ variable).
 epssig: sigP/P max for unpolarized target (sigP/P up to
 epssig doesn’t need standard star when skipdth=True)
 ispol: the Serkowski parameters [P_max, lambda_max, theta_IS]
 to correct IS polarization (P_max ein % and lambda_max in Angstrom). If ispol==None, don’t make the correction of ISP.
 Syntax of out tags: tags1:tags2:tags3, where tags1 is concerning
 to the method to calculate delta_theta correction factor, tags2 is the tag list of the observation of object and tags3 is the tags of the standard star that has been used.
If no standard star was found in some night for one filter, this routine tries to use the standard from another filter to compute the delta theta in missing filter. Case there are no standard star in the other filters also, the routine tries to read a file named as std.link, whose lines content must be [calc night]. It points to a standard from another night night in calcite calc (an example of sintax: [calc night] = [136.1 15out22]). Caso this file was not found, the routine queries directly of what night to use the standard stars.
 Formulas:
 th_eq = fact*th_measured  fact*th^std_measured + th^std_published
 th_eq = fact*th_measured  dth
 dth = fact*th^std_measured  th^std_published
 fact = 1 for observations before 2015 March 1st
 +1 otherwise
 Q_eq = P*cos(2*th_eq*pi/180)
 U_eq = P*sin(2*th_eq*pi/180)
 sigth = 28.65*sig/P

pyhdust.poltools.
graf_qu
(logfile, path2=None, mode=1, thetfile=None, isp=[], odr=True, mcmc=False, nn=[120, 200, 600], thet_ran=[0.0, 180.0], b_ran=[1.0, 1.0], Pb_ran=[0.0, 1.0], Yb_ran=[1.0, 1.0], Vb_ran=[0.0, 1.0], clip=True, sclip=4.5, nmax=5, vfilter=['nostd'], save=False, extens='pdf', limQ=None, limU=None, limJD=None)[source]¶ Plot a QU diagram for the Be star in the logfile .log file (the outfile from polt.genTarget) and fit a line if specified. Propagates error of standard star.
ODR: dotlineMCMC: continuous line
mode=1 plot BVRI graphs in the same figure + U in another; mode=2 plot UBVRI graphs in separated figures.
INPUT
 logfile: Logfile with QU data
 path2: Path to save the output graphs
 mode: 1) Plot a figure to filter U and another for
 BVRI filters; 2) Plot a figure for filter
 thetfile: thet_int.csv file (out from fs.genInt) to
 to plot the lines using the values inside. In this case, mcmc variable don`t matters in the running.
 isp: interstellar polarization to plot direction
 in QU diagram
odr: Run phc.fit_linear to fit a line?
 mcmc: Run fitMCMCline to fit a line?
 nn: fitMCMCline: [n_walkers, n_burnin, n_mcmc]
 thet_ran: fitMCMCline: [thet_min, thet_max]
 b_ran: fitMCMCline: [b_min, b_max]
Pb_ran: fitMCMCline: [Pb_min, Pb_max], with Pb_min >= 0 Yb_ran: fitMCMCline: [Yb_min, Yb_max] Vb_ran: fitMCMCline: [Vb_min, Vb_max], with Vb_min >= 0
clip: phc.fit_linear: apply sigma clipping? sclip: phc.fit_linear: sigma value to clip
 nmax: phc.fit_linear: sigma clipping max number of
 iterations
 vfilter: list of flags to filter (will be marked with
 a ‘x’ symbol and won’t considered in odr fitting). The observations with flag ‘nostd’ never are shown. mcmc fitting uses the filtered observations, except those with ‘nostd’ flag.
save: Save the graphs? If False, just shows
extens: Extension for the graphs
OUTPUT
None if odr and mcmc are False
When mcmc==True: [[[param_u], [sparam_+_u], [sparam__u], n_u, obj_u],
…
[[param_i], [sparam_+_i], [sparam__i], n_i, obj_i]],
where param, sparam_+ and sparam_ are arrays with the five parameters for MCMC fitting (thet, b, Pb, Yb and Vb) and its errors (at right and left), n is the number of points and obj is one graphical dummy object.
When odr==True AND mcmc==False: Idem, but param, sparam_+ and sparam_ are arrays with the two parameters for ODR fitting (a, b) and its errors.
CAUTION:
 The angle returned IS NOT the PA angle obtained from the slop, but the own inclination angle of the fitted line! PA is this angle divided by 2.
 PA angle is indefined by a factor of +90 degrees
 This routine NEVER shows the data with ‘nostd’ tag, independently of vfilter variable!

pyhdust.poltools.
graf_t
(logfile, path2=None, vfilter=['nostd'], save=False, extens='pdf', filt='v')[source]¶ Plot a P_V x t, theta_V x t and P_B/P_I x t graphs for the Be star in the logfile .log file (the outfile from polt.genTarget). Propagates error of standard star.
If ‘nostd’ is in vfilter, no data with ‘nostd’ tag will be displayed, but the others filtered data will be showed with a ‘x’ symbol. If ‘nostd’ is not in vfilter, the data with ‘nostd’ will be displayed normally and the other filtered will be showed with a ‘x’ symbol.

pyhdust.poltools.
grafall
(objdir, filt, nstar=1, n=1, bestouts=[], shortmode=False)[source]¶ Multiple plot modulations for the object inside ‘objdir’ in filter ‘filt’. Return a list ‘sortout’ with the log files sorted by plot number. Allways sortout[0]==’’ and sortout[i] is the plot #i.
 Optional:
 bestouts  list of outs which have smaller error to
 highlight in figures
 shortmode  case True, plot only the groups 16001, 08001,
 08009 (ver .1 and .2) and an eventual 7th group.
 n  is concerning to the nesim 16sequence to use.
 Exemple, if there are 40 WP positions, n=1 is to plot the graphs for positions 116; n=2 is for 1732; n=3 is for 3340.

pyhdust.poltools.
grafpol
(filename, nstar=1, fig=None, ax1=None, ax2=None, save=False, extens='png')[source]¶ Program to plot the best adjust and its residuals of IRAF reduction. ‘filename’ is the path to the .log output file from reduction. nstar is the star number inside out/log files to be plotted.
NEW: Working for *_WP1110….log files! NEW (2): Working for logfiles with more than a single star!
Two working modes:
 If only filename is given: displays the graph or, case save=True, save it into the logfile directory. ‘extens’ parameter changes the extension.
 If filename, one figure and two axes are given: changes these axes, adding plots in them. Doesn’t display, either save at the end, only adds a plot to the axes. Usefull for subplots in same figure.
Authors: Moser and Bednarski Current version: May, 2015

pyhdust.poltools.
lbds
(color, filt, ccdn, airmass=1.3, skiperror=False)[source]¶ Return the lambda_eff in angstrom for star with color index ‘color’, in filter ‘filt’, CCD ‘ccdn’ and airmass ‘airmass’. The default airmass is 1.3 for an average value.
If skiperror==True, tries to use lambda_eff as the value from phc.lbds[] in case of missing information for ‘filt’ or ‘ccd’.
CAUTION: If filt==’u’, the color must be UB
Otherwise, the color must be BV FORMULAS:
 Atm. reddening: redn = 2.5*_np.log10(_np.e)*
 *airmass*(taubtauv)
 Redd. color: color_av = color + redn
 lambda_eff: l_eff = l0 + k1*color_av +
 k2*(color_av**2) +
 k3*(color_av**3)
where tauu, taub, tau are the optical deepth of atmosphere (values used from Kepler de Oliveira et al, Astronomia e Astrofisica).

pyhdust.poltools.
plotMagStar
(tgt, path=None)[source]¶ Function doc
@param PARAM: DESCRIPTION @return RETURN: DESCRIPTION

pyhdust.poltools.
plotfrompollog
(path, star, filters=None, colors=None)[source]¶ Plot default including civil dates

pyhdust.poltools.
propQU
(p, th, sp, sdth, estim='wk')[source]¶ Propagate the delta theta error over the polarization angle, computing the new errors for theta, Q and U.
Input:
 p, th: lists holding the P and theta values.
 sp, sdth: lists holding the P and delta theta errors.
Return lists containing the new errors for theta, Q and U.: sth, sq, su.
Formulas:
sth = sqrt( sth0^2 + sdth^2 ) sq = sqrt( (cos(2*th)*sp)^2 + (p*sin(2*th)*sth)**2 ) su = sqrt( (sin(2*th)*sp)^2 + (p*cos(2*th)*sth)**2 )
Unbias theta error using ‘estim’ estimator:
if p/sp <= K, sth0 = psi otherwise, sth0 = propagated errorwhere K is given by the estimator related to the ‘estim’ variable:
 ‘ml’ : Maximum Likelihood (K=1.41, psi=51.96)
 ‘wk’ : Wardle & Kronberg (K=1.0, psi=51.96)
 ‘’ : None (K=0, psi=51.96)
 ‘mts’: Maier, Tenzer & Santangelo (estimates
 from Bayesian analysis, psi=61.14)

pyhdust.poltools.
queryout
(objdir, obj, f, nstar=1, sigtol=<function <lambda>>)[source]¶ Call chooseout for ‘obj’ at filter ‘f’ (in ‘objdir’), print the graphs and query to the user if he wants to use the selected outs. If not, he musts answer what to use.
 ‘nstar’ == star number inside ‘out’ file (usefull when there are
 more than a single star inside .out)

pyhdust.poltools.
readTests
(tests, tags=None, flag=None)[source]¶ Read boolean list ‘tests’ concerning to dictests dictionary and return a string list with tags and the flag value (‘OK’,’E’,’W’)
‘tags’ and ‘flag’ are optional and are concerning to another tags already assigned and the current flag. The flag returned is the worst flag found between them (e.g., if input ‘flag’ is ‘W’ and tests results on flag ‘OK’, return flag ‘W’; if tests results in flag ‘E’, return ‘E’). Also, if they are given as input, return ‘tags’+tags concerning to ‘tests’ list.

pyhdust.poltools.
readout
(out, nstar=1)[source]¶ Read the *.out file from IRAF reduction and return a float array
Q U SIGMA P THETA SIGMAtheor. APERTURE STAR
 ‘nstar’ == star number inside ‘out’ file (usefull when there are
 more than a single star inside .out)

pyhdust.poltools.
readoutMJD
(out, nstar=1)[source]¶ Read the ‘out’ file from IRAF reduction in a float array (fout), appending the MJD date and the angle of the beams from calcite.
 ‘nstar’ == star number inside ‘out’ file. PS: calcice angle
 is allways evaluated using the first star coordinates.

pyhdust.poltools.
serkowski
(pmax, lmax, wlen, mode, pa=None, law='w82')[source]¶ Mode==1 Receive ISP parameters ‘pmax’, ‘lmax’ e ‘pa’ and return the Stokes QU parameters concerning to the value of Serkowski’s law at wavelenght ‘wlen’.
Mode==2 Receive ISP parameters ‘pmax’, ‘lmax’ and return the P value computed by the Serkowski’s law at wavelenght ‘wlen’.
‘wlen’ and ‘lmax’ must be in Angstrom. ‘wlen’ can be ‘u’/’b’/’v’/’r’/’i’ also.
‘law’ defines what value use to K parameter:
 w82  Wilking (1982)
 K = 1.86*lmax  0.1
 w80  Wilking (1980)
 K = 1.68*lmax  0.002
 serk  Serkowski
 K = 1.15
 Serkowski’s Law:
 P = pmax*np.exp(K*np.log(lmax/wlen)**2)

pyhdust.poltools.
setCCD
(fitsfile)[source]¶ Set CCD name in global variable ‘ccd’.
The CCD name can be: ‘ikon’, ‘ixon’, ‘301’ or ‘ikon14912’ (the last one is the Ikon CCD with Deep Depletion).

pyhdust.poltools.
sintLeff
(ccdn='ixon', step=5.0, save=True, extens='pdf')[source]¶ Sintetizes the response curve, considering the CCD Quantum Efficience (QE) and filter transmitance from OPD and the stellar models from Pickles (1998). Interpolations are made by using cubic splines.
This code DON’T use none curve for sky transmitance!
Creates two data files:
 leff_stars_[ccdn].dat : table with l_eff calculated
 for each star type
 leff_[ccdn].dat : table with the parameters for the
 adjusted cubic function for l_eff(bv) (or l_eff(ub) to the filter U). The Mtype stars were excluded from the adjusts, because the molecular lines were affecting these curves.
 ‘ccdn’: CCD to use in QE curve
 ixon: CON, Frame Transfer ikon: High Sensitive, Frame Transfer
 ikon14912: High Sensitive, Frame Transfer
 301: not avaible
New eventual CCD files must sample the QE, at least, inside the range [2800,11000] Angstrom.
 ‘step’: step, in angstrom, used for the integration
 (Simpson’s method) to calculate lambda_eff. Allowed values are 5, 10, 15, …, 500.
FORMULAS:
The lbd_eff is computed as the fluxweighted mean wavelenght in terms of photons:
lbd_eff = int(lbd*phi(lbd) * d_lbd) / int(phi(lbd) * d_lbd)
 where phi(lbd) = lbd * F_lbd * QE(lbd) * T(lbd)
 F_lbd: stellar flux in erg cm2 s1 A1 QE(lbd): curve for quantum efficience T(lbd): curve for filter transmitance d_lbd: infinitesimal element of wavelength

pyhdust.poltools.
splitData
(night, path_raw='raw', path_red='red')[source]¶ Split the raw files and reduced files for a night.
 Parameters:
 night: path to the night (this directory will be fully preserved) path_raw: directory with the raw data of the nights path_red: directory with the output files of reduction

pyhdust.poltools.
splitData301
(night, path_raw='raw', path_red='red')[source]¶  Split the raw files and reduced files for a night, renaming according
 CCD iXon nomenclature.
 Parameters:
 night: path to the night (this directory will be fully preserved) path_raw: directory with the raw data of the nights path_red: directory with the output files of reduction

pyhdust.poltools.
stdchk
(stdname)[source]¶ Check if the standard star name contains a known name, and return its position in padroes.txt.

pyhdust.poltools.
thtFactor
(MJD)[source]¶ Return the factor for polarization angle ‘theta’. This factor indicates when theta must be taken as 180theta (factor==1) or +theta (factor==+1).
It’s based on engine of IAGPOL polarimeter. If MJD < 57082.5, return 1 (including MJD=1, for no assigned value); otherwise, return +1.
Theta out from polrap is correct for a WP rotating in ‘counterclockwise’ direction, then:
factor = 1 when WP rotating in clockwise factor = +1 when WP rotating in counterclockwise

pyhdust.poltools.
verStdPol
(std, filt, p, sig)[source]¶ Calculate z test for standard ‘std’, filter ‘filt’, comparing the observed polarization ‘p’ + ‘sig’ and the published value.
Return z = abs(ppubp)/sqrt(sigpub^2+sig^2) or 1 if there is no such object or filter.

pyhdust.poltools.
verout
(out, obj, f, nstar=1, verbose=True, delta=3.5)[source]¶ Function to do tests on outfile ‘out’ concerning to star number ‘nstar’, object name ‘obj’ in filter ‘f’. Tests: (1) test if P_pub for standards is compatible with
P_obs value within 10 sigmas (10 because std pol can changes with the time). test if sig < 3*sig_theorical.
 test if there are some standard star (only if
 ‘obj’ is a target star)
 If verbose==True, show warnings in screen
 In objdir==None, outfile is supposed in current dir
Return a boolean list with three components concerning to the tests (1)(3) above + log string. If some test has failed, the concerning value is assigned as “True”; otherwise, “False”.