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
- co-author:
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 burning-in
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 = (1-Pb)/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_i-Yb)**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=['no-std'], 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: dot-line
MCMC: 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 ‘no-std’ never are shown. mcmc fitting uses the filtered observations, except those with ‘no-std’ 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 ‘no-std’ tag, independently of vfilter variable!
- pyhdust.poltools.graf_t(logfile, path2=None, vfilter=['no-std'], 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 ‘no-std’ is in vfilter, no data with ‘no-std’ tag will be displayed, but the others filtered data will be showed with a ‘x’ symbol. If ‘no-std’ is not in vfilter, the data with ‘no-std’ 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 n-esim 16-sequence to use.
Exemple, if there are 40 WP positions, n=1 is to plot the graphs for positions 1-16; n=2 is for 17-32; n=3 is for 33-40.
- 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 U-B
Otherwise, the color must be B-V
- FORMULAS:
- Atm. reddening: redn = 2.5*_np.log10(_np.e)*
*airmass*(taub-tauv)
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 error
where 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 ‘ikon-14912’ (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].dattable with l_eff calculated
for each star type
- leff_[ccdn].dattable with the parameters for the
adjusted cubic function for l_eff(b-v) (or l_eff(u-b) to the filter U). The M-type 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
- ikon-14912: 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 flux-weighted 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 cm-2 s-1 A-1 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 180-theta (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 ‘counter-clockwise’ direction, then:
factor = -1 when WP rotating in clockwise factor = +1 when WP rotating in counter-clockwise
- 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(ppub-p)/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”.