1
2
3 """
4 Tools for trend analysis.
5
6 Author: R. Lombaert
7
8 """
9
10 import os
11 import operator
12 from numpy import array
13 import numpy as np
14 from numpy.random import normal
15 from matplotlib import pyplot as plt
16
17 import cc.path
18 from cc.modeling.objects import Transition
19 from cc.plotting import Plotting2
20 from cc.data import Sed
21
22
23
24 -def makeDiagnosticPlot(sg,molec,scaling=[],escaling=[],combine_water=0,\
25 edists=[],pfn_path=''):
26
27 '''
28 Make a diagnostic plot for a series of stars, where Imb for transitions of
29 a given molecule are compared with their upper level energies.
30
31 Line strengths are always scaled with distance squared. Additional scaling
32 can be requested.
33
34 Three plots are made: One with scaling versus distance and two with scaling
35 versus distance and the CO line strength of the J=15-14 and J=30-29 lines
36 respectively. The comparison with CO line strengths is not scaled except
37 with distance.
38
39 @param sg: The stellar models, in which the transitions have been matched
40 with integrated line strengths.
41 @type sg: list[Star()]
42 @param molec: The molecule for which this is done, shorthand notation.
43 @type molec: str
44
45 @keyword scaling: Scale the line strengths also with a keyword given here.
46 'MDOT_GAS', etc. Assumes a Star() object knows the
47 keyword.
48
49 (default:[])
50 @type scaling: list[string]
51 @keyword combine_water: Combine both ortho and para water in a single plot
52
53 (default: 0)
54 @type combine_water: bool
55 @keyword edists: Include errors for distance estimates of stars here.
56
57 (default: [])
58 @type edists: list[float]
59 @keyword escaling: Include relative errors of extra scaling parameters
60 as lists. len is len of scaling, len of an element
61 is len of sg
62
63 (default: [])
64 @type escaling: list[list]
65 @keyword pfn_path: Output folder for diagnostic plots. Default if to be
66 stored locally.
67
68 (default: '')
69 @type pfn_path: str
70
71 '''
72
73 if '1H1H16O' not in molec: combine_water = 0
74 if len(scaling) != len(escaling):
75 print 'No errors on the scaling factors taken into account.'
76 escaling = []
77 if len(edists) != len(sg):
78 print 'No errors on the distance taken into account.'
79 edists = [0]*len(sg)
80 allints = []
81 allenergies = []
82 allerrors = []
83 allwaves = []
84 allints_co15 = []
85 allints_co30 = []
86 allerrs_co15 = []
87 allerrs_co30 = []
88
89
90 co = sg[0].getMolecule('12C16O')
91
92 co1514 = Transition.Transition(molecule=co,telescope='PACS',jup=15,\
93 jlow=14,path_gastronoom='codeJun2013')
94 co3029 = Transition.Transition(molecule=co,telescope='PACS',jup=30,\
95 jlow=29,path_gastronoom='codeJun2013')
96 trl_co15 = Transition.getTransFromStarGrid(sg,co1514,'sample')
97 trl_co30 = Transition.getTransFromStarGrid(sg,co3029,'sample')
98 ls_co15,errs_co15 = Transition.getLineStrengths(trl_co15,'dint')
99 ls_co30,errs_co30 = Transition.getLineStrengths(trl_co30,'dint')
100
101 for istar,(s,intco15,intco30,errco15,errco30) in \
102 enumerate(zip(sg,ls_co15,ls_co30,errs_co15,errs_co30)):
103 sn = s['STAR_NAME']
104 allints.append([])
105 allenergies.append([])
106 allerrors.append([])
107 allwaves.append([])
108 allints_co15.append([])
109 allints_co30.append([])
110 allerrs_co15.append([])
111 allerrs_co30.append([])
112 if combine_water: \
113 trans = s.getTransitions('1H1H16O') + s.getTransitions('p1H1H16O')
114 else:
115 trans = s.getTransitions(molec)
116
117 ls_trans,errls_trans = Transition.getLineStrengths(trans,'dint')
118 for t,lst,elst in zip(trans,ls_trans,errls_trans):
119 if not np.isnan(lst):
120 allints_co15[-1].append(abs(lst/intco15))
121 allints_co30[-1].append(abs(lst/intco30))
122 allerrs_co15[-1].append(np.sqrt(sum([elst**2]+[errco15**2]))\
123 *abs(lst/intco15))
124 allerrs_co30[-1].append(np.sqrt(sum([elst**2]+[errco30**2]))\
125 *abs(lst/intco30))
126
127 tint = abs(lst*s['DISTANCE']**2/100.**2 )
128 for par in scaling:
129 tint *= 1./s[par]
130 totalerr = np.sqrt(sum([elst**2]+\
131 [2*edists[istar]**2]+\
132 [esca[istar]**2 for esca in escaling]))
133 allints[-1].append(tint)
134 allenergies[-1].append(t.getEnergyUpper())
135 allerrors[-1].append(totalerr*tint)
136 allwaves[-1].append(t.wavelength*10**4)
137 else:
138 for tlist in [allints_co15,allints_co30,allerrs_co15,\
139 allerrs_co30,allints,allerrors]:
140 tlist[-1].append(float('nan'))
141 allenergies[-1].append(float(t.getEnergyUpper()))
142 allwaves[-1].append(t.wavelength*10**4)
143 isort = np.argsort(allenergies[-1])
144 for tlist in [allints_co15,allints_co30,allerrs_co15,allwaves,\
145 allerrs_co30,allints,allerrors,allenergies]:
146 tlist[-1] = array(tlist[-1])[isort]
147
148 isort = np.argsort([s['MDOT_GAS'] for s in sg])
149 for tlist in [allints_co15,allints_co30,allerrs_co15,allwaves,\
150 allerrs_co30,allints,allerrors,allenergies,sg]:
151 tlist = array(tlist)[isort]
152
153 linestyles = ['o-','o-','o-','o-','o-','o-','o-',\
154 '--x','--x','--x','--x','--x','--x','--x',\
155 '-.s','-.s','-.s','-.s','-.s','-.s','-.s']
156 colors = ['r','b','k','g','m','y','c']
157 line_types = [ls + col for ls,col in zip(linestyles,3*colors)]
158 line_types = line_types[:len(sg)]
159
160 xmin = 0
161 xmax = len(allenergies[0])+1
162 plot_title = ', '.join(['%i: %.1f cm$^{-1}$ - %.1f $\mu$m'\
163 %(i+1,t.getEnergyUpper(),t.wavelength*10**4)
164 for i,t in enumerate(sg[0].getTransitions(molec))])
165 if scaling:
166 plot_title += 'Extra scaling: ' + ', '.join([sca.replace('_','\_')
167 for sca in scaling])
168 x = []
169 y = []
170 yerr = []
171 y_co15 = []
172 ye15 = []
173 y_co30 = []
174 ye30 = []
175 for istar in range(len(sg)):
176 x.append([])
177 y.append([])
178 yerr.append([])
179 y_co15.append([])
180 y_co30.append([])
181 ye15.append([])
182 ye30.append([])
183 for i,(iint,eint,iint_co15,iint_co30,eint_co15,eint_co30) in \
184 enumerate(zip(allints[istar],allerrors[istar],\
185 allints_co15[istar],allints_co30[istar],\
186 allerrs_co15[istar],allerrs_co30[istar])):
187 if not np.isnan(iint):
188 x[-1].append(i+1)
189 y[-1].append(iint)
190 yerr[-1].append(eint)
191 if not np.isnan(iint_co15):
192 y_co15[-1].append(iint_co15)
193 ye15[-1].append(eint_co15)
194 if not np.isnan(iint_co30):
195 y_co30[-1].append(iint_co30)
196 ye30[-1].append(eint_co30)
197 extension = 'pdf'
198
199 path = os.path.join(pfn_path,'ints_vs_eul_%s'%molec)
200 for par in scaling:
201 if par == scaling[0]:
202 path += '_extrascaling'
203 path += '_%s'%par
204 print Plotting2.plotCols(filename=path,\
205 x=x,y=y,extension=extension,\
206 yerr=yerr,\
207 xmin=xmin,xmax=xmax,plot_title=plot_title,\
208 yaxis='$I_\mathrm{int}$ (W/m$^2$)',\
209 xaxis='Index Energy Upper Level',\
210 line_types=line_types,\
211 keytags=['%.1e -- %s'%(s['MDOT_GAS'],s['STAR_NAME']) for s in sg],\
212 key_location=(0.87,0.01),ylogscale=1,\
213 linewidth=3,fontsize_key=26,fontsize_title=20)
214
215 if molec != '12C16O' and not scaling:
216 path = os.path.join(pfn_path,'ints_co15_vs_eul_%s'%molec)
217 print Plotting2.plotCols(filename=path,\
218 x=x,y=y_co15,extension=extension,\
219 yerr=ye15,\
220 xmin=xmin,xmax=xmax,plot_title=plot_title,\
221 yaxis='$I_\mathrm{int}/I_\mathrm{CO 15-14}$ (W/m$^2$)',\
222 xaxis='Index Energy Upper Level',\
223 line_types=line_types,\
224 keytags=['%.1e -- %s'%(s['MDOT_GAS'],s['STAR_NAME']) for s in sg],\
225 key_location=(0.87,0.01),ylogscale=1,\
226 linewidth=3,fontsize_key=26,fontsize_title=20)
227
228 path = os.path.join(pfn_path,'ints_co30_vs_eul_%s'%molec)
229 print Plotting2.plotCols(filename=path,\
230 x=x,y=y_co30,extension=extension,\
231 yerr=ye30,\
232 xmin=xmin,xmax=xmax,plot_title=plot_title,\
233 yaxis='$I_\mathrm{int}/I_\mathrm{CO 30-29}$ (W/m$^2$)',\
234 xaxis='Index Energy Upper Level',\
235 line_types=line_types,\
236 keytags=['%.1e -- %s'%(s['MDOT_GAS'],s['STAR_NAME']) for s in sg],\
237 key_location=(0.87,0.01),ylogscale=1,\
238 linewidth=3,fontsize_key=26,fontsize_title=20)
239
240
241
242 -def makeParamPlot(sg,xpar,ypar,expar=[],eypar=[],xratios=[],yratios=[],\
243 emdot=[],exparlog=0,eyparlog=0,edists=[],mode='dint',\
244 n_data=0,extra_mpar=[],extra_dpar=[],cfg='',pfn_path='',\
245 add_linear_fit=0,alf_xmin=None,alf_xmax=None,seds=[],\
246 deredden=0,**kwargs):
247
248 '''
249 Make a diagnostic plot of either measured line strengths or intrinsic
250 parameters versus measured line strengths or intrinsic parameters.
251
252 Ratios are possible for line strengths. Not for intrinsic parameters.
253
254 Requires preparatory work done for the Pacs() and the Star() objects.
255
256 @param sg: The stellar models, in which the transitions have been matched
257 with integrated line strengths. If both models and data are
258 combined, the data Star() objects are assumed to be listed
259 first.
260 @type sg: list[Star()]
261 @param xpar: The parameter on the x-axis. Can be either a string (Star()
262 keyword), or an index (of the transition in the first object
263 in the sg list) for line strengths, or a float giving the
264 wavelength of the continuum point. When looking at line
265 strengths in a combo mode (cint or
266 ctmb) this means it is the index in the transition list of the
267 data objects rather than the model objects. Transitions in
268 objects other than the first 1 can have different indices.
269 Note the essential difference between floats and integers!
270 @type xpar: string/int
271 @param ypar: The parameter on the y-axis. Can be either a string (Star()
272 keyword), or an index (of the transition in the first object
273 in the sg list) for line strengths, or a float giving the
274 wavelength of the continuum point. When looking at line
275 strengths inn a combo mode (cint or
276 ctmb) this means it is the index in the transition list of the
277 data objects rather than the model objects. Transitions in
278 objects other than the first 1 can have different indices.
279 Note the essential difference between floats and integers!
280 @type ypar: string/int
281
282 @keyword xratios: If xpar is a line strength or a continuum point, multiple
283 ratios can be requested to be plotted in succession.
284 Therefore, this gives the indices (if int, refers to the
285 1st Star() object in sg) or 'mdot' (if ratio wrt Mdot)
286 or float (in case of a continuum wavelength point) for the
287 x-axis ratio.
288
289 (default: [])
290 @type xratios: list[int/str]
291 @keyword yratios: If ypar is a line strength, multiple ratios can be
292 requested to be plotted in succession.
293 Therefore, this gives the indices (if int, refers to the
294 1st Star() object in sg) or 'mdot' (if ratio wrt Mdot)
295 or float (in case of a continuum wavelength point) for the
296 y-axis ratio
297
298 (default: [])
299 @type yratios: list[int/str]
300 @keyword emdot: Include errors for the x/yratio quantity if it is mdot. Not
301 used for error estimation on mdot as a parameter! The mdot
302 errors are given in log scale.
303
304 (default: [])
305 @type emdot: list[float]
306 @keyword expar: The error on the x-parameter if it is a Star() key and if
307 mode is cint or dint. Number of entries in array is equal
308 to the number of data Star() objects.
309
310 (default: [])
311 @type expar: array
312 @keyword eypar: The error on the y-parameter if it is a Star() key and if
313 mode is cint or dint. Number of entries in array is equal
314 to the number of data Star() objects.
315
316 (default: [])
317 @type eypar: array
318 @keyword exparlog: The xpar error is given in logscale. Only relevant for
319 the d and c modes.
320
321 (default: 0)
322 @type exparlog: bool
323 @keyword eyparlog: The ypar error is given in logscale. Only relevant for
324 the d and c modes.
325
326 (default: 0)
327 @type eyparlog: bool
328 @keyword edists: Include errors for distance estimates of stars here. These
329 distances are only used to rescale line strengths if they
330 are not in a ratio.
331
332 (default: [])
333 @type edists: list[float]
334 @keyword mode: The mode in which line strengths are selected, ie either
335 from data or from models. Either 'dint', 'mint', 'mtmb' or
336 'dtmb' values. A combination of both is possible by setting
337 this key to 'cint' or 'ctmb'. Then the extra keyword
338 'n_data' is required, which indicates how many Star()
339 objects are associated with data. The method assumes they
340 are the first objects in the list of Star() objects. In case
341 only continuum wavelengths are requested (a combination is
342 possible!), only the first letter really matters.
343
344 (default: 'dint')
345 @type mode: str
346 @keyword n_data: The number of data Star() objects, assuming they are the
347 first in the star_grid. Only required if mode == 'combo'.
348 This number, if given, must be equal to the number of seds,
349 if given.
350
351 (default: 0)
352 @type n_data: int
353 @keyword extra_mpar: If extra conditional parameters are requested for
354 models, the
355 plot is colour coded based on them. For instance,
356 Star() object keywords can serve as conditionals.
357 Note that these are only applied when mode == mtmb,
358 mint, cint or ctmb.
359
360 (default: [])
361 @type extra_mpar: list[string]
362 @keyword extra_dpar: If extra conditional parameters are requested for data,
363 the plot
364 is colour coded based on them. For instance, Star()
365 object keywords can serve as conditionals. Note that
366 these are only applied when mode == dtmb, dint, cint
367 or ctmb.
368
369 (default: [])
370 @type extra_dpar: list[string]
371 @keyword seds: The SEDs of the data objects. Only used when xpar or ypar is
372 a float (and thus continuum points are required).
373 The number of SEDs given must
374 be equal to n_data, or the number of Star() objects if
375 mode[0] == 'd'. An error is thrown otherwise.
376
377 (default: [])
378 @type seds: list[Sed()]
379 @keyword deredden: Deredden the SEDs before plotting, in case of continuum
380 flux points. This is never done in case only one data
381 object is given and reddening is requested in models to
382 avoid double correction.
383
384 (default: 0)
385 @type deredden: bool
386 @keyword cfg: config filename read as a dictionary, can replace any keyword
387 given to plotCols. Can also be a dictionary itself, in which
388 case no file is read and the kwargs are updated with the
389 content of cfg
390
391 (default: '')
392 @type cfg: string/dict
393 @keyword add_linear_fit: Add a linear fit to the figures. The fit is done
394 through corrSG method, of which extra arguments
395 can be given in kwargs. (xpar_co, ypar_co)
396 Only works in dint or cint mode if xratios or
397 yratios has len less than 2.
398
399 (default: 0)
400 @type add_linear_fit: bool
401 @keyword pfn_path: Output folder for diagnostic plots. Default if to be
402 stored locally.
403
404 (default: '')
405 @type pfn_path: str
406 @keyword alf_xmin: The minimum x value for the linear fit plot, if
407 requested. (This is not the cut off value for the
408 fitting routine itself!) Has to be given if a linear fit
409 is requested.
410
411 (default: None)
412 @type alf_xmin: float
413 @keyword alf_xmax: The maximum x value for the linear fit plot, if
414 requested. (This is not the cut off value for the
415 fitting routine itself!) Has to be given if a linear fit
416 is requested.
417
418 (default: None)
419 @type alf_xmax: float
420
421 @keyword **kwargs: extra keywords needed for the linear fit, if requested.
422 @type **kwargs: dict
423
424 @return: The filename of the produced plot is returned.
425 @rtype: str
426
427 '''
428
429 x_titles = dict([('MDOT_GAS',r'$\log$ $\left[\dot{M}_\mathrm{g}\ (\mathrm{M}_\odot/\mathrm{yr})\right]$'),\
430 ('MDOT_DUST',r'$\log$ $\left[\dot{M}_\mathrm{d}\ (\mathrm{M}_\odot/\mathrm{yr})\right]$'),\
431 ('VEL_INFINITY_GAS',r'$v_{\infty\mathrm{,g}}$ ($\mathrm{km} \mathrm{s}^{-1}$)'),\
432 ('SHELLMASS',r'$\log$ $\left[\dot{M}_\mathrm{g}/v_{\infty\mathrm{,g}}\ (\mathrm{M}_\odot\ \mathrm{yr}^{-1}\ \mathrm{km}^{-1}\ \mathrm{s})\right]$'),\
433 ('SHELLDENS',r'$\log$ $\left[\bar{\rho}\ (\mathrm{g}\ \mathrm{cm}^{-3})\right]$'),\
434 ('SHELLCOLDENS',r'$\log$ $\left[\bar{m}\ (\mathrm{g}\ \mathrm{cm}^{-2})\right]$'),\
435 ('SHELLDENS2',r'$\sqrt{\bar{\rho}^2 R_\star}$ (g/cm$^{5/2}$)'),\
436 ('L_STAR','$L_\star$ (L$_\odot$)'),\
437 ('P_STAR',r'$\log$ $\left[P\ (\mathrm{days})\right]$'),\
438 ('T_STAR','$T_\star$ (K)'),\
439 ('R_STAR','$R_\star$ (Rsun)'),\
440 ('Q_STAR','$Q_\star$ (days)'),\
441 ('R_INNER_GAS','$R_\mathrm{i,g}$ (R$_\star$)'),\
442 ('AH2O_RATE',r'$\log$ $\left[A_{\mathrm{H}_2\mathrm{O}}/A_{\mathrm{H}_2} \times \dot{M}_\mathrm{g}\ (\mathrm{M}_\odot/\mathrm{yr})\right]$'),\
443 ('F_H2O',r'$\log$ $\left[A_{\mathrm{H}_2\mathrm{O}}/A_{\mathrm{H}_2}\right]$'),\
444 ])
445 pfn_parts = dict([('MDOT_GAS','mg'),\
446 ('MDOT_DUST',r'md'),\
447 ('STARTYPE','startype'),\
448 ('A_SICB','asicb'),\
449 ('A_AMCSPH','aamcsph'),\
450 ('VEL_INFINITY_GAS','vg'),\
451 ('SHELLMASS','shellmass'),\
452 ('SHELLDENS','dens'),\
453 ('SHELLCOLDENS','coldens'),\
454 ('SHELLDENS2','dens3-2'),\
455 ('L_STAR','lstar'),\
456 ('R_STAR','rstar'),\
457 ('P_STAR','period'),\
458 ('Q_STAR','qstar'),\
459 ('T_STAR','tstar'),\
460 ('P_TYPE','ptype'),\
461 ('F_H2O','ah2o'),\
462 ('DUST_TO_GAS_CHANGE_ML_SP','d2g'),\
463 ('TEMPERATURE_EPSILON_GAS','eps1'),\
464 ('TEMPERATURE_EPSILON2_GAS','eps2'),\
465 ('TEMPERATURE_EPSILON3_GAS','eps3'),\
466 ('RADIUS_EPSILON2_GAS','rt12'),\
467 ('RADIUS_EPSILON3_GAS','rt23'),\
468 ('R_INNER_GAS','rig'),\
469 ('MDOT_CLASS','mdotgrad'),\
470 ('SCD_CLASS','scdgrad'),\
471 ('SHELLMASS_CLASS','shellmassclass'),\
472 ('ABUN_O','abuno'),\
473 ('L_CLASS','lclass'),\
474 ('T_CLASS','tclass'),\
475 ('VG_CLASS','vgclass'),\
476 ('AH2O_RATE','ah2orate'),\
477 ('F_CONT_TYPE','fconttype'),\
478 ('DRIFT_TYPE','drifttype'),\
479 ('ENHANCE_ABUNDANCE_FACTOR_H2O','h2oabunfac'),\
480 ('ABUNDANCE_FILENAME_H2O','h2oabunfile')])
481 keynames = dict([('MDOT_GAS','$\dot{M}_\mathrm{g}$'),\
482 ('MDOT_DUST',r'$\dot{M}_\mathrm{d}$'),\
483 ('A_SICB','A(SICB)'),\
484 ('STARTYPE','StarType'),\
485 ('A_AMCSPH','A(AMCSPH)'),\
486 ('VEL_INFINITY_GAS',r'$v_{\infty\mathrm{,g}}$'),\
487 ('SHELLMASS',r'$\dot{M}_\mathrm{g}/v_{\infty\mathrm{,g}}$'),\
488 ('SHELLDENS',r'$\bar{\rho}$'),\
489 ('SHELLCOLDENS',r'$\bar{m}$'),\
490 ('SHELLDENS2',r'$\sqrt{\bar{\rho}^2 R_\star}}$'),\
491 ('L_STAR','$L_\star$'),\
492 ('P_STAR','$P$'),\
493 ('Q_STAR','$Q_\star$'),\
494 ('T_STAR','$T_\star$'),\
495 ('P_TYPE','Var.~Type'),\
496 ('F_H2O','$n_{\mathrm{H}_2\mathrm{O}}/n_{\mathrm{H}_2}$'),\
497 ('DUST_TO_GAS_CHANGE_ML_SP','$\psi$'),\
498 ('TEMPERATURE_EPSILON_GAS',r'$\epsilon$'),\
499 ('TEMPERATURE_EPSILON2_GAS',r'$\epsilon_2$'),\
500 ('TEMPERATURE_EPSILON3_GAS',r'$\epsilon_3$'),\
501 ('RADIUS_EPSILON2_GAS','$R_\mathrm{T, 12}$'),\
502 ('RADIUS_EPSILON3_GAS','$R_\mathrm{T, 23}$'),\
503 ('R_INNER_GAS','$R_\mathrm{i,g}$'),\
504 ('MDOT_CLASS',''),\
505 ('SCD_CLASS',''),\
506 ('SHELLMASS_CLASS',''),\
507 ('ABUN_O','$n_{\mathrm{O}}/n_{\mathrm{H}_\mathrm{tot}}$'),\
508 ('L_CLASS',''),\
509 ('T_CLASS',''),\
510 ('VG_CLASS',''),\
511 ('DRIFT_TYPE',''),\
512 ('AH2O_RATE',r'$A_{\mathrm{H}_2\mathrm{O}}/A_{\mathrm{H}_2} \times \dot{M}_\mathrm{g}$'),\
513 ('F_CONT_TYPE','Type F$_\mathrm{cont}$'),\
514 ('ENHANCE_ABUNDANCE_FACTOR_H2O','h2oAbunFac'),\
515 ('ABUNDANCE_FILENAME_H2O','h2oAbunFile')])
516 keyunits = dict([('MDOT_GAS','$\mathrm{M}_\odot\ \mathrm{yr}^{-1}$'),\
517 ('MDOT_DUST','$\mathrm{M}_\odot\ \mathrm{yr}^{-1}$'),\
518 ('STARTYPE',''),\
519 ('A_SICB',''),\
520 ('A_AMCSPH',''),\
521 ('VEL_INFINITY_GAS','$\mathrm{km\;s}^{-1}$'),\
522 ('SHELLMASS','$\mathrm{M}_\odot\ \mathrm{yr}^{-1}\ \mathrm{km}^{-1}\ \mathrm{s}$'),\
523 ('SHELLDENS','$\mathrm{g\;cm}^{-3}$'),\
524 ('SHELLCOLDENS','$\mathrm{g\;cm}^{-2}$'),\
525 ('SHELLDENS2','$\mathrm{g\;cm}^{5/2}$'),\
526 ('L_STAR','$\mathrm{L}_\odot$'),\
527 ('P_STAR','$\mathrm{days}$'),\
528 ('Q_STAR','$\mathrm{days}$'),\
529 ('T_STAR','$\mathrm{K}$'),\
530 ('P_TYPE',''),\
531 ('F_H2O',''),\
532 ('DUST_TO_GAS_CHANGE_ML_SP',''),\
533 ('TEMPERATURE_EPSILON_GAS',''),\
534 ('TEMPERATURE_EPSILON2_GAS',''),\
535 ('TEMPERATURE_EPSILON3_GAS',''),\
536 ('RADIUS_EPSILON2_GAS','$\mathrm{R}_\star$'),\
537 ('RADIUS_EPSILON3_GAS','$\mathrm{R}_\star$'),\
538 ('R_INNER_GAS','$\mathrm{R}_\star$'),\
539 ('MDOT_CLASS',''),\
540 ('SCD_CLASS',''),\
541 ('SHELLMASS_CLASS',''),\
542 ('ABUN_O',''),\
543 ('T_CLASS',''),\
544 ('VG_CLASS',''),\
545 ('DRIFT_TYPE','Drift'),\
546 ('L_CLASS',''),\
547 ('AH2O_RATE','$\mathrm{M}_\odot\ \mathrm{yr}^{-1}$'),\
548 ('F_CONT_TYPE',''),\
549 ('ENHANCE_ABUNDANCE_FACTOR_H2O',''),\
550 ('ABUNDANCE_FILENAME_H2O','')])
551 makeints = dict([('MDOT_GAS',0),\
552 ('MDOT_DUST',0),\
553 ('STARTYPE',0),\
554 ('A_SICB',0),\
555 ('A_AMCSPH',0),\
556 ('VEL_INFINITY_GAS',1),\
557 ('SHELLMASS',0),\
558 ('SHELLDENS',0),\
559 ('SHELLCOLDENS',0),\
560 ('SHELLDENS2',0),\
561 ('L_STAR',1),\
562 ('P_STAR',1),\
563 ('Q_STAR',0),\
564 ('T_STAR',1),\
565 ('P_TYPE',0),\
566 ('F_H2O',0),\
567 ('DUST_TO_GAS_CHANGE_ML_SP',0),\
568 ('TEMPERATURE_EPSILON_GAS',0),\
569 ('TEMPERATURE_EPSILON2_GAS',0),\
570 ('TEMPERATURE_EPSILON3_GAS',0),\
571 ('RADIUS_EPSILON2_GAS',0),\
572 ('RADIUS_EPSILON3_GAS',0),\
573 ('R_INNER_GAS',0),\
574 ('MDOT_CLASS',0),\
575 ('SCD_CLASS',0),\
576 ('SHELLMASS_CLASS',0),\
577 ('ABUN_O',0),\
578 ('L_CLASS',0),\
579 ('T_CLASS',0),\
580 ('VG_CLASS',0),\
581 ('DRIFT_TYPE',0),\
582 ('AH2O_RATE',0),\
583 ('F_CONT_TYPE',0),\
584 ('ENHANCE_ABUNDANCE_FACTOR_H2O',0),\
585 ('ABUNDANCE_FILENAME_H2O',0)])
586
587 for k in extra_mpar+extra_dpar:
588 if k not in pfn_parts.keys():
589 pfn_parts[k] = k.lower().replace('_','')
590 keynames[k] = k.replace('_','\_')
591 keyunits[k] = ''
592 makeints[k] = 0
593
594 edists,emdot = array(edists), array(emdot)
595 n_data = int(n_data)
596 expar, eypar = array(expar), array(eypar)
597 if isinstance(extra_mpar,str):
598 extra_mpar = [extra_mpar]
599 if isinstance(extra_dpar,str):
600 extra_dpar = [extra_dpar]
601
602
603 if isinstance(xpar,str):
604 xratios = []
605 if isinstance(ypar,str):
606 yratios = []
607 ratios = [xratios,yratios]
608 sg_dists = array([s['DISTANCE'] for s in sg])
609 sg_mdot = array([s['MDOT_GAS'] for s in sg])
610
611
612
613 pars = [xpar,ypar]
614 epars = [expar,eypar]
615 eparlogs = [exparlog,eyparlog]
616
617
618
619
620 if isinstance(xpar,float) or isinstance(ypar,float) \
621 or True in [isinstance(i,float) for i in xratios+yratios]:
622 add_linear_fit = 0
623 if mode[0] == 'd' and not seds:
624 raise IOError('No SEDs given for data objects.')
625 elif mode[0] == 'c' and n_data != len(seds):
626 raise IOError('Number of SEDs not equal to number of data objects')
627
628
629
630 if mode[0] == 'm':
631
632 add_linear_fit = 0
633 expar = array([])
634 eypar = array([])
635 n_data = 0
636 seds = []
637 extra_dpar = []
638 extra_par = extra_mpar
639 current_par = 'm'
640 elif mode[0] == 'd':
641 n_data = len(sg)
642 extra_mpar = []
643 if mode[0] != 'm':
644 extra_par = extra_dpar
645 current_par = 'd'
646 for istar,s in enumerate(sg):
647 if n_data and istar == n_data:
648 extra_par = extra_mpar
649 current_par = 'm'
650 s['EC'] = (current_par,\
651 tuple([s[par]
652 for par in extra_par
653 if not s[par] is None]))
654 ecl = sorted(list(set([s['EC'] for s in sg])))
655 ecl_num = []
656 for ec in ecl:
657 isg = array([s['EC'] == ec for s in sg])
658 isgd = isg[:n_data]
659 isgm = isg[n_data:]
660 nsgd = len(isg[:n_data][isgd])
661 nsgm = len(isg[n_data:][isgm])
662 ecl_num.append((ec,isg,isgd,isgm,nsgd,nsgm))
663
664 if len(xratios) > 1 or len(yratios) > 1:
665
666 add_linear_fit = 0
667 if add_linear_fit:
668 add_linear_fit = dict()
669 xratio, yratio = None, None
670 if xratios:
671 xratio = xratios[0]
672 if yratios:
673 yratio = yratios[0]
674 results = corrSG(sg=sg[:n_data],xpar=xpar,ypar=ypar,expar=expar,\
675 eypar=eypar,xratio=xratio,yratio=yratio,edist=edists,\
676 show=0,eyratio=eyratio,exratio=exratio,**kwargs)
677 fitcoef = results['results']
678 xgrid = []
679 ygrid = []
680 this_x = array([alf_xmin,alf_xmax])
681 add_linear_fit['xmean'] = this_x
682 add_linear_fit['ymean'] = results['intercept']+results['slope']*this_x
683 for n in range(0, fitcoef.shape[0],4):
684 this_y = fitcoef[n,1] + fitcoef[n,0] * this_x
685 xgrid.append(this_x)
686 ygrid.append(this_y)
687 add_linear_fit['xgrid'] = xgrid
688 add_linear_fit['ygrid'] = ygrid
689
690
691
692
693
694 if 'mdot' in xratios:
695 xratios[xratios.index('mdot')] = 'xmdot'
696 if 'mdot' in yratios:
697 yratios[yratios.index('mdot')] = 'ymdot'
698
699
700
701 ls_ratios = dict()
702 els_ratios = dict()
703 for i in set(yratios+xratios):
704
705 if i == 'xmdot' or i == 'ymdot': continue
706 elif isinstance(i,float):
707
708
709 dists = [s['DISTANCE'] for s in sg[:n_data]] if deredden else []
710 cflux,eflux = Sed.getCFlux(wav=i,seds=seds,star_grid=sg[n_data:],
711 deredden=dists)
712 ls_ratios[i] = cflux
713 els_ratios[i] = eflux
714 else:
715 ratsample = sg[0]['GAS_LINES'][i]
716 rattrans = Transition.getTransFromStarGrid(sg,ratsample,'sample')
717 ls,els = Transition.getLineStrengths(rattrans,mode,n_data=n_data)
718 ls_ratios[i] = ls
719 els_ratios[i] = els
720
721
722
723
724 x, xerr, xblend = dict([('x',[])]), dict([('x',[])]), dict([('x',[])])
725 for k in xratios:
726 x[k], xerr[k], xblend[k] = [], [], []
727 y, yerr, yblend = dict([('y',[])]), dict([('y',[])]), dict([('y',[])])
728 for k in yratios:
729 y[k], yerr[k], yblend[k] = [], [], []
730
731 blends = [xblend,yblend]
732 xy = [x,y]
733 errs = [xerr,yerr]
734 axes = ['x','y']
735
736
737
738
739
740
741
742 sample, seltrans, allint, allerr = [] , [], [], []
743 for par,epar,eparlog,blend,xyi,err,axisstr in zip(pars,epars,eparlogs,\
744 blends,xy,errs,axes):
745 if isinstance(par,str):
746 sg_par = array([s[par] for s in sg])
747 sample.append(None)
748 seltrans.append(None)
749 allint.append(None)
750 allerr.append(None)
751 for (ec,isg,isgd,isgm,nsgd,nsgm) in ecl_num:
752
753 blend[axisstr].append(np.zeros(nsgd) != 0)
754 xyi[axisstr].append(np.log10(sg_par[isg]))
755
756 if epar.size and eparlog:
757 err[axisstr].append(np.concatenate([epar[isgd],np.zeros(nsgm)]))
758
759 elif epar.size and len(epar.shape) == 1:
760 ll = np.concatenate([-np.log10(1-epar[isgd]),np.zeros(nsgm)])
761 ul = np.concatenate([np.log10(1+epar[isgd]),np.zeros(nsgm)])
762 err[axisstr].append([ll,ul])
763
764 elif epar.size and epar.shape[0] == 2:
765 ll = np.concatenate([-np.log10(1-epar[0][isgd]),np.zeros(nsgm)])
766 ul = np.concatenate([np.log10(1+epar[1][isgd]),np.zeros(nsgm)])
767 err[axisstr].append([ll,ul])
768 elif isinstance(par,float):
769
770
771 dists = [d for d in sg_dists[:n_data]] if deredden else []
772 all_iflux, all_ieflux = Sed.getCFlux(wav=par,seds=seds,\
773 star_grid=sg[n_data:],\
774 deredden=dists)
775 sample.append(None)
776 seltrans.append(None)
777 allint.append(all_iflux)
778 allerr.append(all_ieflux)
779 else:
780
781
782
783
784
785 isample = sg[0]['GAS_LINES'][par]
786 iseltrans = Transition.getTransFromStarGrid(sg,isample,'sample')
787 iallint,iallerr = Transition.getLineStrengths(iseltrans,mode,\
788 n_data=n_data)
789 sample.append(isample)
790 seltrans.append(iseltrans)
791 allint.append(iallint)
792 allerr.append(iallerr)
793
794
795
796
797
798 for iratios,iallint,iallerr,axisstr in zip(ratios,allint,allerr,axes):
799 if axisstr+'mdot' in iratios:
800 line1 = abs(iallint[:n_data])*sg_dists[:n_data]**2/100.**2
801 line1_err = line1*np.sqrt(iallerr[:n_data]**2+4*edists**2)
802 line2 = np.log10(sg_mdot[:n_data])
803 line2_err = emdot
804 mratios = guessRatio(line1,line1_err,line2,line2_err,\
805 line2_log=1,positive=1,n_fit=10000)
806 emrat = array([np.std(mratios[:,istar])/np.mean(mratios[:,istar])
807 for istar in range(len(line1))])
808 ls_ratios[axisstr+'mdot'] = sg_mdot
809 els_ratios[axisstr+'mdot'] = emrat
810
811
812
813 for (ec,isg,isgd,isgm,nsgd,nsgm) in ecl_num:
814
815 for par,blend,xyi,err,axisstr,iallint,iallerr,iratios \
816 in zip(pars,blends,xy,errs,axes,allint,allerr,ratios):
817 if not isinstance(par,str):
818
819 blend[axisstr].append(iallint[isgd] < 0)
820
821 irat1 = abs(iallint[isg])
822 xyi[axisstr].append(np.log10(irat1*sg_dists[isg]**2/100.**2))
823
824
825
826
827
828
829 if mode[0] != 'm':
830 etot = np.sqrt(iallerr[isgd]**2+4*edists[isgd]**2)
831 ll = np.concatenate([-np.log10(1-etot),np.zeros(nsgm)])
832 ul = np.concatenate([np.log10(1+etot),np.zeros(nsgm)])
833 err[axisstr].append([ll,ul])
834
835 for k in iratios:
836 if k == axisstr+'mdot':
837
838 blend[k].append(blend[axisstr][-1])
839 xyi[k].append(xyi[axisstr][-1]-np.log10(ls_ratios[k][isg]))
840 if mode[0] != 'm':
841 etot = els_ratios[k][isgd]
842 else:
843
844
845 blend[k].append((blend[axisstr][-1])+(ls_ratios[k][isgd]<0))
846 xyi[k].append(np.log10(irat1/abs(ls_ratios[k][isg])))
847 if mode[0] != 'm':
848 etot = np.sqrt(iallerr[isgd]**2+els_ratios[k][isgd]**2)
849 if mode[0] != 'm':
850 ll = np.concatenate([-np.log10(1-etot),np.zeros(nsgm)])
851 ul = np.concatenate([np.log10(1+etot),np.zeros(nsgm)])
852 err[k].append([ll,ul])
853
854
855 cfg_dict = Plotting2.readCfg(cfg)
856 if cfg_dict.has_key('pfn_path'):
857 pfn_path = cfg_dict['pfn_path']
858 extra_ppars = dict()
859
860 pt = ''
861 if isinstance(ypar,int):
862 pt += '%s: E$_\mathrm{ul,y}$ = %.1f - %.2f'\
863 %(str(sample[1]),sample[1].getEnergyUpper(),\
864 sample[1].wavelength*10**4)
865 if isinstance(xpar,int):
866 pt += 'VS %s: E$_\mathrm{ul,x}$ = %.1f - %.2f'\
867 %(str(sample[0]),sample[0].getEnergyUpper(),\
868 sample[0].wavelength*10**4)
869 extra_ppars['plot_title'] = pt
870 extra_ppars['fontsize_title'] = 20
871 extra_ppars['figsize'] = (8*np.sqrt(2),8)
872 extra_ppars['extension'] = '.pdf'
873 extra_ppars['fontsize_key'] = 14
874 extra_ppars['linewidth'] = 2
875
876
877
878 if not ecl in [[('m',()),('d',())],[('m',())],[('d',())]]:
879 keytags = []
880 for curr_par,ec in ecl:
881 this_par = curr_par == 'm' and extra_mpar or extra_dpar
882 k = []
883 for par,v in zip(this_par,ec):
884 if 'CLASS' in par:
885 kstr = v[1]
886 elif par == 'P_TYPE':
887 kstr = '$\mathrm{%s}$'%v
888 elif par == 'SHELLCOLDENS':
889 kstr = '%s = $%.2f$ %s'%(keynames[par],v,keyunits[par])
890 else:
891 kstr = '%s = $%s$ %s'\
892 %(keynames[par],\
893 makeints[par] and str(int(v)) or str(v),\
894 keyunits[par])
895 k.append(kstr)
896 keytags.append(', '.join(k))
897 extra_ppars['key_location'] = 'best'
898 mlinestyles = ['-x','-x','-x','-x','-x','-x','-x',\
899 '--s','--s','--s','--s','--s','--s','--s',\
900 '-.+','-.+','-.+','-.+','-.+','-.+','-.+',\
901 '--p','--p','--p','--p','--p','--p','--p',\
902 'o-','o-','o-','o-','o-','o-','o-']
903 dlinestyles = ['o','o','o','o','o','o','o',\
904 'x','x','x','x','x','x','x',\
905 's','s','s','s','s','s','s']
906 colors = ['r','b','g','k','m','y','c']
907 dline_types = [ls + col for ls,col in zip(dlinestyles,3*colors)]
908 colors.reverse()
909 mline_types = [ls + col for ls,col in zip(mlinestyles,5*colors)]
910 if mode[0] == 'm':
911 line_types = mline_types[:len(ecl)]
912 zorder = range(len(ecl))
913 elif mode[0] == 'd':
914 line_types = dline_types[:len(ecl)]
915 zorder = range(len(ecl))
916 else:
917 d_ecl = len([ec for ec in ecl if ec[0] == 'd'])
918 m_ecl = len([ec for ec in ecl if ec[0] == 'm'])
919 line_types = dline_types[:d_ecl] \
920 + mline_types[:m_ecl]
921 zorder = range(10,10+d_ecl) + range(-m_ecl,0)
922 markersize = [6]*len(keytags)
923 pfn_ecl = '_'+'_'.join([pfn_parts[ec] for ec in extra_dpar+extra_mpar])
924 if pfn_ecl == '_': pfn_ecl = ''
925
926
927
928
929 if xratios:
930 del x['x']
931 if yratios:
932 del y['y']
933
934
935 for xk in x.keys():
936
937 if isinstance(xk,int):
938 xratsample = sg[0]['GAS_LINES'][xk]
939
940
941
942
943 if not xerr[xk]:
944 extra_ppars['xmin'] = min([min(xi) for xi in x[xk]])-0.2
945 extra_ppars['xmax'] = max([max(xi) for xi in x[xk]])+0.2
946
947
948 if isinstance(xpar,str):
949 extra_ppars['xaxis'] = x_titles[xpar]
950 pfn_xtag = pfn_parts[xpar]
951 pfn_xrat = ''
952 elif isinstance(xpar,float):
953 s1 = r'F_\mathrm{%.1f\ \mu m}'%xpar
954 if xk == 'xmdot':
955 extra_ppars['xaxis'] = r'$\log$ $\left[%s/\dot{M}_\mathrm{'%s1+\
956 r'g}\ (\mathrm{W}/\mathrm{m}^2/\mathrm{Hz}\ '+\
957 r'\mathrm{yr}/\mathrm{M}_\odot)\right]$'
958 elif xk == 'x':
959 extra_ppars['xaxis'] = r'$\log$ $\left[%s\ (\mathrm{W}/\mathrm{m}^2/\mathrm{Hz})\right]$'%s1
960 elif isinstance(xk,int):
961 s2 = xratsample.makeAxisLabel()
962 extra_ppars['xaxis'] = r'$\log$ $\left[%s/%s\ (\mathrm{Hz}^{-1})\right]$'%(s1,s2)
963 else:
964 s2 = r'F_\mathrm{%.1f\ \mu m}'%xk
965 extra_ppars['xaxis'] = r'$\log$ $\left[%s/%s\right]$'%(s1,s2)
966 pfn_xtag = 'f%.1fmic'%xpar
967 if xk =='xmdot':
968 pfn_xrat = '_mdot'
969 elif xk == 'x':
970 pfn_xrat = ''
971 elif isinstance(xk,int):
972 ms = yratsample.molecule.molecule_short
973 pfn_xrat = '_%s%i'%(ms,xk)
974 else:
975 pfn_xrat = '_f%.1fmic'%xk
976 else:
977
978 s1 = sample[0].makeAxisLabel()
979 if xk == 'xmdot':
980 extra_ppars['xaxis'] = r'$\log$ $\left[%s/\dot{M}_\mathrm{'%s1+\
981 r'g}\ (\mathrm{W}/\mathrm{m}^2\ \mathrm{yr}/\mathrm{M}_\odot)\right]$'
982 elif xk == 'x':
983 extra_ppars['xaxis'] = r'$\log$ $\left[%s\ (\mathrm{W}/\mathrm{m}^2)\right]$'%s1
984 elif isinstance(xk,float):
985 s2 = r'F_\mathrm{%.1f\ \mu m}'%xk
986 extra_ppars['xaxis'] = r'$\log$ $\left[%s/%s'%(s1,s2)+\
987 r'\ (\mathrm{Hz})\right]$'
988 else:
989 iml = sample[0].molecule.molecule != xratsample.molecule.molecule
990 s1 = sample[0].makeAxisLabel(iml)
991 s2 = xratsample.makeAxisLabel(iml)
992 extra_ppars['xaxis'] = r'$\log$ $\left[%s/%s\right]$'%(s1,s2)
993
994 pfn_xtag = '%s_eul_%i_wl_%.1f'\
995 %(sample[0].molecule.molecule,\
996 int(sample[0].getEnergyUpper()),\
997 float(sample[0].wavelength*10**4))
998 if xk =='xmdot':
999 pfn_xrat = '_mdot'
1000 elif xk == 'x':
1001 pfn_xrat = ''
1002 elif isinstance(xk,float):
1003 pfn_xrat = '_f%.1fmic'%xk
1004 else:
1005 ms = xratsample.molecule.molecule_short
1006 pfn_xrat = '_%s%i'%(ms,xk)
1007
1008
1009 for yk in y.keys():
1010 if isinstance(yk,int):
1011 yratsample = sg[0]['GAS_LINES'][yk]
1012
1013
1014
1015
1016 if not yerr[yk]:
1017 extra_ppars['ymin'] = min([min(yi) for yi in y[yk]])-0.2
1018 extra_ppars['ymax'] = max([max(yi) for yi in y[yk]])+0.2
1019
1020 if isinstance(ypar,str):
1021 extra_ppars['yaxis'] = x_titles[ypar]
1022 pfn_ytag = pfn_parts[ypar]
1023 pfn_yrat = ''
1024 elif isinstance(ypar,float):
1025 s1 = r'F_\mathrm{%.1f\ \mu m}'%ypar
1026 if yk == 'ymdot':
1027 extra_ppars['yaxis'] = r'$\log$ $\left[%s/\dot{M}_\mathrm{'%s1+\
1028 r'g}\ (\mathrm{W}/\mathrm{m}^2/\mathrm{Hz}\ '+\
1029 r'\mathrm{yr}/\mathrm{M}_\odot)\right]$'
1030 elif yk == 'y':
1031 extra_ppars['yaxis'] = r'$\log$ $\left[%s\ (\mathrm{W}/\mathrm{m}^2/\mathrm{Hz})\right]$'%s1
1032 elif isinstance(yk,int):
1033 s2 = yratsample.makeAxisLabel()
1034 extra_ppars['yaxis'] = r'$\log$ $\left[%s/%s\ (\mathrm{Hz}^{-1})\right]$'%(s1,s2)
1035 else:
1036 s2 = r'F_\mathrm{%.1f\ \mu m}'%yk
1037 extra_ppars['yaxis'] = r'$\log$ $\left[%s/%s\right]$'%(s1,s2)
1038 pfn_ytag = 'f%.1fmic'%ypar
1039 if yk =='ymdot':
1040 pfn_yrat = '_mdot'
1041 elif yk == 'y':
1042 pfn_yrat = ''
1043 elif isinstance(yk,int):
1044 ms = yratsample.molecule.molecule_short
1045 pfn_yrat = '_%s%i'%(ms,yk)
1046 else:
1047 pfn_yrat = '_f%.1fmic'%yk
1048 else:
1049 s1 = sample[1].makeAxisLabel()
1050 if yk == 'ymdot':
1051 extra_ppars['yaxis'] = r'$\log$ $\left[%s/\dot{M}_\mathrm{'%s1+\
1052 r'g}\ (\mathrm{W}/\mathrm{m}^2\ \mathrm{yr}/\mathrm{M}_\odot)\right]$'
1053 elif yk == 'y':
1054 extra_ppars['yaxis'] = r'$\log$ $\left[%s\ (\mathrm{W}/\mathrm{m}^2)\right]$'%s1
1055 elif isinstance(yk,float):
1056 s2 = r'F_\mathrm{%.1f\ \mu m}'%yk
1057 extra_ppars['yaxis'] = r'$\log$ $\left[%s/%s'%(s1,s2)+\
1058 r'\ (\mathrm{Hz})\right]$'
1059 else:
1060 iml = sample[1].molecule.molecule != yratsample.molecule.molecule
1061 s1 = sample[1].makeAxisLabel(iml)
1062 s2 = yratsample.makeAxisLabel(iml)
1063 extra_ppars['yaxis'] = r'$\log$ $\left[%s/%s\right]$'%(s1,s2)
1064
1065 pfn_ytag = '%s_eul_%i_wl_%.1f'\
1066 %(sample[1].molecule.molecule,\
1067 int(sample[1].getEnergyUpper()),\
1068 float(sample[1].wavelength*10**4))
1069 if yk =='ymdot':
1070 pfn_yrat = '_mdot'
1071 elif yk == 'y':
1072 pfn_yrat = ''
1073 elif isinstance(yk,float):
1074 pfn_yrat = '_f%.1fmic'%yk
1075 else:
1076 ms = yratsample.molecule.molecule_short
1077 pfn_yrat = '_%s%i'%(ms,yk)
1078
1079
1080 xb, yb = [[]], [[]]
1081 for blend1,blend2,xi,yi in zip(xblend[xk],yblend[yk],x[xk],y[yk]):
1082 blended = blend1 + blend2
1083 xb[-1].extend(xi[blended])
1084 yb[-1].extend(yi[blended])
1085
1086 if xb[-1]:
1087 extra_ppars['keytags'] = keytags + ['$\mathrm{Blended}$']
1088 extra_ppars['line_types'] = line_types + ['xk']
1089 extra_ppars['markersize'] = markersize + [14]
1090 extra_ppars['zorder'] = zorder + [max(zorder)+1]
1091 else:
1092 xb, yb, = [], []
1093 extra_ppars['keytags'] = keytags
1094 extra_ppars['line_types'] = line_types
1095 extra_ppars['markersize'] = markersize
1096 extra_ppars['zorder'] = zorder
1097
1098 if add_linear_fit:
1099 extra_ppars['keytags'] = extra_ppars['keytags'] + ['Mean Linear fit']
1100 extra_ppars['line_types'] = extra_ppars['line_types'] + ['-g'] + ['-k']*len(add_linear_fit['xgrid'])
1101 extra_ppars['markersize'] = extra_ppars['markersize'] + [4] + [4]*len(add_linear_fit['xgrid'])
1102 extra_ppars['zorder'] = extra_ppars['zorder'] + [min(zorder)-1] + [min(zorder)-2]*len(add_linear_fit['xgrid'])
1103 extra_ppars['alpha'] = [1]*len(x[xk])+[1]*len(xb)+[1]+[0.002]*len(add_linear_fit['xgrid'])
1104 xb.append(add_linear_fit['xmean'])
1105 yb.append(add_linear_fit['ymean'])
1106 xb.extend(add_linear_fit['xgrid'])
1107 yb.extend(add_linear_fit['ygrid'])
1108
1109
1110
1111
1112 extra_ppars.update(cfg_dict)
1113 pfn = os.path.join(pfn_path,'%s_%s%s_vs_%s%s%s'\
1114 %(mode,pfn_ytag,pfn_yrat,pfn_xtag,\
1115 pfn_xrat,pfn_ecl))
1116
1117 extra_ppars['filename'] = pfn
1118 ff = Plotting2.plotCols(x=xb and x[xk]+xb or x[xk],\
1119 y=yb and y[yk]+yb or y[yk],\
1120 yerr=yb and yerr[yk]+[None]*len(yb) or yerr[yk],\
1121 xerr=xb and xerr[xk]+[None]*len(xb) or xerr[xk],\
1122 **extra_ppars)
1123 print ff
1124
1125 if n_data > 0:
1126 print 'Stars plotted (in order of x):'
1127 for xi,yi,(ec,isg,isgd,isgm,nsgd,nsgm) \
1128 in zip(x[xk],y[yk],ecl_num):
1129 if ec[0] == 'm': continue
1130 ifin = np.isfinite(xi) * np.isfinite(yi)
1131 isort = np.argsort(xi[ifin])
1132 sgsort = array(sg)[isgd][ifin][isort]
1133 k = ', '.join(['%s = %s%s%s'
1134 %(keynames[con],\
1135 makeints[con] and str(int(v)) or str(v),\
1136 keyunits[con] and ' ' or '',\
1137 keyunits[con])
1138 for con,v in zip(extra_dpar,ec[1])])
1139 print k, ': %s'%', '.join([s['STAR_NAME'] for s in sgsort])
1140
1141 return ff
1142
1143
1144 -def guessRatio(line1,line1_err,line2,line2_err,line1_log=0,line2_log=0,\
1145 n_fit=10000,positive=0):
1146
1147 '''
1148 Guess a ratio of given values with error bars a given number of times.
1149
1150 For both components of the ratio, values are drawn from a Gaussian
1151 distribution around the given value, with the error as sigma.
1152
1153 Can be used for error analysis, or for estimating errors on ratios in case
1154 the errors on the separate components are difficult to propagate properly.
1155
1156 The option to guess a value within error bars in log space is possible. The
1157 resulting value is then converted back to linear space, after which the
1158 ratio is taken.
1159
1160 Negative values can occur, due to the Gaussian nature of the guesses. Keep
1161 this in mind when taking the log of output values. If you do not want
1162 negative values, this van be requested via the positive keyword.
1163
1164 A guess of the ratio, and a standard deviation, can be calculated by taking
1165 the mean and std of the columns in the ouput array.
1166
1167 @param line1: Values of the first parameter on the y-axis
1168 @type line1: array
1169 @param line1_err: Uncertainties on line1, assuming they are in a normal
1170 distribution (1-sigma)
1171 @type line1_err: array
1172 @param line2: Values of the second parameter on the y-axis. Can be an empty
1173 array in case you want to fit a correlation between par and
1174 line1 without any ratio involved. Pass an empty array if you
1175 simply want to randomize a single array, instead of a ratio.
1176 @type line2: array
1177 @param line2_err: Uncertainties on line2, assuming they are in a normal
1178 distribution (1-sigma)
1179 @type line2_err: array
1180 @keyword line1_log: If line 1 is in log scale. In the ratio, 10**line1 is
1181 then taken.
1182
1183 (default: 0)
1184 @type line1_log: bool
1185 @keyword line2_log: if line 2 is in log scale. In the ratio, 10**line1 is
1186 then taken.
1187
1188 (default: 0)
1189 @type line2_log: bool
1190 @keyword n_fit: The number of times the correlation is fitted.
1191
1192 (default: 10000)
1193 @type n_fit: int
1194 @keyword positive: In some cases, you may want to disallow negative values,
1195 eg when you take the log of the results. This switch
1196 allows you to exclude negative values from the output.
1197 Use this with caution! In some case, this will severely
1198 affect the Gaussian distribution.
1199
1200 (default: 0)
1201 @type positive: bool
1202
1203 @return: The n_fit guesses of the requested ratio.
1204 @rtype: array((n_fit,len(line1)))
1205
1206 '''
1207
1208 n_fit = int(n_fit)
1209 line1_log, line2_log = bool(line1_log), bool(line2_log)
1210 line1, line1_err = array(line1), array(line1_err)
1211 line2, line2_err = array(line2), array(line2_err)
1212 yarr = np.empty((n_fit,len(line1)))
1213 if positive:
1214 for n in range(n_fit):
1215 while True:
1216 guess1 = normal(line1, line1_err)
1217 if line1_log:
1218 guess1 = 10**guess1
1219 break
1220 elif False not in (guess1[np.isfinite(guess1)] > 0):
1221 break
1222
1223 while True:
1224 if line2.size != 0:
1225 guess2 = normal(line2, line2_err)
1226 else:
1227 guess2 = 1
1228 break
1229 if line2_log:
1230 guess2 = 10**guess2
1231 break
1232 elif False not in (guess2[np.isfinite(guess2)] > 0):
1233 break
1234 yarr[n] = guess1/guess2
1235 else:
1236 for n in range(n_fit):
1237 guess1 = normal(line1, line1_err)
1238 if line1_log:
1239 guess1 = 10**guess1
1240 if line2.size != 0:
1241 guess2 = normal(line2, line2_err)
1242 else:
1243 guess2 = 1
1244 if line2_log and line2.size != 0:
1245 guess2 = 10**guess2
1246 yarr[n] = guess1/guess2
1247 return yarr
1248
1249
1250 -def fitCorrPolyLog(par1,par1_err,par2,par2_err,line1,line1_err,line2,line2_err,\
1251 par1_log=0,par2_log=0,line1_log=0,line2_log=0,n_fit=10000,\
1252 poly_degree=1,show=0,fn_plt='',x_for_yratio=0,\
1253 y_for_xratio=0):
1254
1255 '''
1256 Fit a polynomial to a data set.
1257
1258 The data set can consist of straight up values or of ratios on both the x
1259 and y axis (in log space).
1260
1261 Takes into account errors in both dimensions.
1262
1263 Can be used for e.g. error estimation on a correlation.
1264
1265 @param par1: Values of the first parameter on x-axis.
1266 @type par1: array
1267 @param par1_err: Uncertainties on par, assuming they are in a normal
1268 distribution (1-sigma)
1269 @type par1_err: array
1270 @param par2: Values of the second parameter on the x-axis. Can be an empty
1271 array in case you don't want a ratio on the x-axis.
1272 @type par2: array
1273 @param par2_err: Uncertainties on par2, assuming they are in a normal
1274 distribution (1-sigma)
1275 @type par2_err: array
1276 @param line1: Values of the first parameter on the y-axis
1277 @type line1: array
1278 @param line1_err: Uncertainties on line1, assuming they are in a normal
1279 distribution (1-sigma)
1280 @type line1_err: array
1281 @param line2: Values of the second parameter on the y-axis. Can be an empty
1282 array in case you don't want a ratio on the x-axis.
1283 @type line2: array
1284 @param line2_err: Uncertainties on line2, assuming they are in a normal
1285 distribution (1-sigma)
1286 @type line2_err: array
1287 @keyword par1_log: If par is in log scale. If not, the log will be taken of
1288 par, since this method fits log log correlations.
1289
1290 (default: 0)
1291 @type par1_log: bool
1292 @keyword par2_log: If par2 is in log scale. If not, the log will be taken
1293 of par2, since this method fits log log correlations.
1294
1295 (default: 0)
1296 @type par2_log: bool
1297 @keyword line1_log: If line 1 is in log scale. In the ratio, 10**line1 is
1298 then taken.
1299
1300 (default: 0)
1301 @type line1_log: bool
1302 @keyword line2_log: if line 2 is in log scale. In the ratio, 10**line1 is
1303 then taken.
1304
1305 (default: 0)
1306 @type line2_log: bool
1307 @keyword n_fit: The number of times the correlation is fitted.
1308
1309 (default: 10000)
1310 @type n_fit: int
1311 @keyword poly_degree: The degree of the polynomial that is fitted.
1312
1313 (default: 1)
1314 @type poly_degree: int
1315 @keyword show: Show a plot with the results. If cfg is given, the plot is
1316 adapted, including the filename.
1317
1318 (default: 0)
1319 @type show: bool
1320 @keyword fn_plt: The filename of the plot, in case show is True, and a
1321 saved plot is requested.
1322
1323 (default: '')
1324 @type fn_plt: str
1325 @keyword x_for_yratio: Use the par grid as the second component in the y
1326 ratio. This can be useful for instance if the ratio
1327 has Mdot as numerator, while Mdot is also on the x
1328 axis. In this case, you want to use the same random
1329 value for the same point on both x and y.
1330
1331 (default: 0)
1332 @type x_for_yratio: bool
1333 @keyword y_for_xratio: Use the line1 grid as the second component in the x
1334 ratio. This can be useful for instance if the ratio
1335 has Mdot as numerator, while Mdot is also on the y
1336 axis. In this case, you want to use the same random
1337 value for the same point on both x and y.
1338
1339 (default: 0)
1340 @type y_for_xratio: bool
1341
1342 @return: The fit results are returned for all n_fit fitted functions. The
1343 parameters are the output of np.polyfit and the amount depends on
1344 the polynomial degree.
1345 @rtype: array
1346
1347 '''
1348
1349 poly_degree = int(poly_degree)
1350
1351 fitcoef = np.empty((n_fit, poly_degree+1))
1352 if y_for_xratio:
1353 xarr = guessRatio(par1,par1_err,[],[],line1_log=par1_log,n_fit=n_fit,\
1354 positive=1)
1355 y1 = guessRatio(line1,line1_err,[],[],line1_log=line1_log,n_fit=n_fit,\
1356 positive=1)
1357 else:
1358 xarr = guessRatio(par1,par1_err,par2,par2_err,line1_log=par1_log,\
1359 line2_log=par2_log,n_fit=n_fit,positive=1)
1360
1361 if x_for_yratio:
1362 yarr = guessRatio(line1,line1_err,[],[],line1_log,\
1363 n_fit=n_fit,positive=1)
1364 x1 = guessRatio(par1,par1_err,[],[],line1_log=par1_log,\
1365 n_fit=n_fit,positive=1)
1366 else:
1367 yarr = guessRatio(line1,line1_err,line2,line2_err,line1_log,line2_log,\
1368 n_fit=n_fit,positive=1)
1369
1370 for n,x,y in zip(range(n_fit),xarr,yarr):
1371
1372
1373
1374
1375
1376
1377 xl = np.log10(x)
1378 yl = np.log10(y)
1379 if x_for_yratio:
1380 yl = yl - np.log10(x1[n])
1381 if y_for_xratio:
1382 xl = xl - np.log10(y1[n])
1383 fitcoef[n] = np.polyfit(xl, yl, poly_degree)
1384
1385 if show and poly_degree == 1:
1386
1387 plt.figure(1)
1388 plt.clf()
1389
1390
1391 plt.subplot(221)
1392
1393 x1 = par1
1394 if par2.size != 0:
1395 x2 = par2
1396 else:
1397 x2 = 1
1398 if par1_log:
1399 x1 = 10**x1
1400 if par2_log and par2.size != 0:
1401 x2 = 10**x2
1402 x = np.log10(x1/x2)
1403
1404 y1 = line1
1405 if line2.size != 0:
1406 y2 = line2
1407 else:
1408 y2 = 1
1409 if line1_log:
1410 y1 = 10**y1
1411 if line2_log and line2.size != 0:
1412 y2 = 10**y2
1413 y = np.log10(y1/y2)
1414
1415 plt.scatter(x, y, color='blue', marker='o')
1416
1417 x_grid = np.linspace(1.05*x.min(),0.95*x.max(),100)
1418 for n in range(0, fitcoef.shape[0], 2):
1419 y_grid = fitcoef[n,1] + fitcoef[n,0] * x_grid
1420 plt.plot(x_grid, y_grid, color="red", alpha = 0.006)
1421
1422
1423 plt.xlabel("log(X)")
1424 plt.ylabel("log(Y)")
1425
1426 plt.subplot(222)
1427 plt.hexbin(fitcoef[:,0], fitcoef[:,1], bins=40)
1428 plt.xlabel("Slope")
1429 plt.ylabel("Intercept")
1430
1431 plt.subplot(223)
1432 plt.hist(fitcoef[:,0], bins=40)
1433 plt.xlabel("Slope")
1434 plt.ylabel("N")
1435
1436 plt.subplot(224)
1437 plt.hist(fitcoef[:,1], bins=40)
1438 plt.xlabel("Intercept")
1439 plt.ylabel("N")
1440
1441 if not fn_plt:
1442 plt.show()
1443 else:
1444 plt.savefig(fn_plt)
1445
1446 return fitcoef
1447
1448
1450
1451 '''
1452 Tool for selecting data from a star_grid. Mainly used by corrStarGrid() to
1453 define input arrays for the correlation study.
1454
1455 @param sg: The grid of Star() objects.
1456 @type sg: list[Star()]
1457 @param par: The requested parameter. If a string, a Star() keyword is used
1458 and if MDOT_GAS then the error bars are set as log10(3)/3. If
1459 an integer, it is the index of a Transition() in the Star()
1460 and the line strength and -error on it- of the line is taken.
1461 @type par: str/int
1462 @param epar: The 1-sigma error bar of the parameter. Only relevant if par
1463 is a string Star() keyword that is not MDOT_GAS.
1464 @type epar: array
1465 @keyword par_co: Define cutoff values here. Always given as an array of
1466 size 2. If a lower and/or upper boundary is not needed,
1467 it is set as None. In case of MDOT_GAS==xpar, these values
1468 are converted to log scale. Can be set as None if no
1469 cutoff is needed. Only relevant if par is a string.
1470
1471 (default: None)
1472 @type par_co: array
1473
1474 @keyword edist: Give the relative error of the distance here. Used to
1475 estimate an uncertainty on the rescaled line strengths
1476 according to distance (down to 100 pc). Not relevant when
1477 par is a string. An empty array implies no scaling.
1478
1479 (default: [])
1480 @type edist: array
1481
1482 @return: The values for the requested parameter, as well as the
1483 uncertainty, the cutoff values and the log request, ie all
1484 keywords needed for the fitCorrPolyLog method.
1485 @rtype: (arr,arr,arr,bool)
1486
1487 '''
1488
1489 edist = array(edist)
1490 if isinstance(par,str):
1491 vals = array([s[par] for s in sg])
1492 if not par_co is None:
1493 if par_co[0] is None:
1494 par_co[0] = min(vals)
1495 if par_co[1] is None:
1496 par_co[1] = max(vals)
1497 par_co = array(par_co,dtype=np.float)
1498 if par == 'MDOT_GAS':
1499 evals = np.ones(len(sg))*np.log10(3.)/3.
1500 vals = np.log10(vals)
1501 if not par_co is None: par_co = np.log10(par_co)
1502 vals_log = 1
1503 else:
1504 evals = epar
1505 vals_log = 0
1506 else:
1507 sample = sg[0]['GAS_LINES'][par]
1508 trans = Transition.getTransFromStarGrid(sg,sample,'sample')
1509 vals,evals = Transition.getLineStrengths(trans,mode='dint')
1510
1511
1512
1513 vals = abs(vals)
1514 vals_log = 0
1515 if edist.size:
1516 dists = array([s['DISTANCE'] for s in sg])
1517 vals = vals*dists**2/100.**2
1518 evals = vals*np.sqrt(4*edist**2+evals**2)
1519 else:
1520 evals = vals*evals
1521
1522 finvals = vals[np.isfinite(vals)]
1523 par_co = array([min(finvals),max(finvals)])
1524 return (vals,evals,par_co,vals_log)
1525
1526
1527
1528 -def corrSG(sg,xpar,ypar,expar=[],eypar=[],xratio=None,yratio=None,\
1529 eyratio=[],exratio=[],edist=[],xpar_co=(None,None),\
1530 ypar_co=(None,None),**kwargs):
1531
1532 '''
1533 A method focused on finding correlations between parameters and/or data
1534 of multiple Star() objects.
1535
1536 @param sg: The grid of Star() objects.
1537 @type sg: list[Star()]
1538 @param xpar: The parameter on x-axis. If a string, a Star() keyword is used
1539 and if MDOT_GAS then the error bars are set as log10(3)/3. If
1540 an integer, it is the index of a Transition() in the Star()
1541 and the line strength and -error on it- of the line is taken.
1542 If the same as yratio, the same guesses are used for both.
1543 @type xpar: str/int
1544 @param ypar: The parameter on y-axis. If a string, a Star() keyword is used
1545 and if MDOT_GAS then the error bars are set as log10(3)/3. If
1546 an integer, it is the index of a Transition() in the Star()
1547 and the line strength and -error on it- of the line is taken.
1548 If the same as xratio, the same guesses are used for both.
1549 @type ypar: int/str
1550
1551 @keyword expar: The 1-sigma error bar of the parameter on the xaxis. Only
1552 relevant if xpar is a string Star() keyword that is not
1553 MDOT_GAS.
1554
1555 (default: [])
1556 @type expar: array
1557 @keyword eypar: The 1-sigma error bar of the parameter on the yaxis. Only
1558 relevant if ypar is a string Star() keyword that is not
1559 MDOT_GAS.
1560
1561 (default: [])
1562 @type eypar: array
1563 @keyword xratio: If a ratio on the x-axis is requested, this is the second
1564 component. Input syntax same as xpar.
1565
1566 (default: None)
1567 @type xratio: int/str
1568 @keyword exratio: The relative error bars on the xratio parameter are given
1569 here. Only relevant if xratio is a string Star() keyword
1570 that is not MDOT_GAS.
1571
1572 (default: [])
1573 @type exratio: array
1574 @keyword yratio: If a ratio on the y-axis is requested, this is the second
1575 component. Input syntax same as ypar.
1576
1577 (default: None)
1578 @type yratio: int/str
1579 @keyword eyratio: The relative error bars on the yratio parameter are given
1580 here. Only relevant if yratio is a string Star() keyword
1581 that is not MDOT_GAS.
1582
1583 (default: [])
1584 @type eyratio: array
1585 @keyword edist: Give the relative error of the distance here. Used to
1586 estimate an uncertainty on the rescaled line strengths
1587 according to distance (down to 100 pc). Not relevant when a
1588 line ratio is requested.
1589
1590 (default: [])
1591 @type edist: array
1592 @keyword xpar_co: Define cutoff values here. Always given as an array of
1593 size 2. If a lower and/or upper boundary is not needed,
1594 it set as None. In case of MDOT_GAS==xpar, these values
1595 are converted to log scale.
1596
1597 (default: array(None,None))
1598 @type xpar_co: array
1599 @keyword ypar_co: Define cutoff values here. Always given as an array of
1600 size 2. If a lower and/or upper boundary is not needed,
1601 it set as None. In case of MDOT_GAS==ypar, these values
1602 are converted to log scale.
1603
1604 (default: array(None,None))
1605 @type ypar_co: array
1606 @keyword kwargs: Extra keywords relevant for fitLinearCorr or
1607 fitCorrPolyLog
1608 @type kwargs: dict
1609
1610 @return: The resulting fit parameters are returned. NYI if poly_degree!=1.
1611 @rtype: dict()
1612
1613 '''
1614
1615
1616
1617
1618 expar, eypar, edist = array(expar), array(eypar), array(edist)
1619 exratio, eyratio = array(exratio), array(eyratio)
1620 xpar_co, ypar_co = array(xpar_co), array(ypar_co)
1621 ep = dict()
1622
1623
1624 if xratio is None or isinstance(xratio,str):
1625 xv, exv, xpar_co, ep['par1_log'] = selectDataSG(sg,xpar,expar,xpar_co,\
1626 edist)
1627 else:
1628 xv, exv, xpar_co, ep['par1_log'] = selectDataSG(sg,xpar,expar,xpar_co)
1629
1630
1631 if xratio is None or xratio == ypar:
1632 if not xratio is None: ep['y_for_xratio'] = 1
1633 xrat, exrat, ep['par2_log'] = array([]),array([]),0
1634 else:
1635 xrat, exrat, dummy, ep['par2_log'] = selectDataSG(sg,xratio,exratio)
1636
1637
1638 if yratio is None or isinstance(yratio,str):
1639 yv, eyv, ypar_co, ep['line1_log'] = selectDataSG(sg,ypar,eypar,\
1640 ypar_co,edist)
1641 else:
1642 yv, eyv, ypar_co, ep['line1_log'] = selectDataSG(sg,ypar,eypar,ypar_co)
1643
1644
1645 if yratio is None or yratio == xpar:
1646 if not yratio is None: ep['x_for_yratio'] = 1
1647 yrat, eyrat, ep['line2_log'] = array([]),array([]),0
1648 else:
1649 yrat, eyrat, dummy, ep['line2_log'] = selectDataSG(sg,yratio,eyratio)
1650
1651
1652 isort = np.argsort(xv)
1653 xv, exv = xv[isort], exv[isort]
1654 yv, eyv = yv[isort], eyv[isort]
1655 if xrat.size:
1656 xrat, exrat = xrat[isort], exrat[isort]
1657 if yrat.size:
1658 yrat, eyrat = yrat[isort], eyrat[isort]
1659
1660
1661 if xrat.size:
1662 xselfinite = np.isfinite(xv/xrat)
1663 else:
1664 xselfinite = np.isfinite(xv)
1665 if yrat.size:
1666 yselfinite = np.isfinite(yv/yrat)
1667 else:
1668 yselfinite = np.isfinite(yv)
1669 selfinite = yselfinite * xselfinite
1670 if xrat.size:
1671 xrat, exrat = xrat[selfinite], exrat[selfinite]
1672 if yrat.size:
1673 yrat, eyrat = yrat[selfinite], eyrat[selfinite]
1674 xv, exv = xv[selfinite], exv[selfinite]
1675 yv, eyv = yv[selfinite], eyv[selfinite]
1676
1677
1678 xbools = (xv>=xpar_co[0]) * (xv<=xpar_co[1])
1679 ybools = (yv>=ypar_co[0]) * (yv<=ypar_co[1])
1680 bools = xbools*ybools
1681 xv, exv = xv[bools], exv[bools]
1682 yv, eyv = yv[bools], eyv[bools]
1683 if xrat.size:
1684 xrat, exrat = xrat[bools], exrat[bools]
1685 if yrat.size:
1686 yrat, eyrat = yrat[bools], eyrat[bools]
1687 kwargs.update(ep)
1688
1689 allcoef = fitCorrPolyLog(par1=xv,par1_err=exv,par2=xrat,par2_err=exrat,\
1690 line1=yv,line1_err=eyv,line2=yrat,\
1691 line2_err=eyrat,**kwargs)
1692 results = dict()
1693 if kwargs.get('poly_degree',1) == 1:
1694 results['n_points'] = len(xv)
1695 results['slope'] = allcoef[:,0].mean()
1696 results['eslope'] = allcoef[:,0].std()
1697 results['intercept'] = allcoef[:,1].mean()
1698 results['eintercept'] = allcoef[:,1].std()
1699
1700 for x, x_err, name in [(results['slope'],results['eslope'],"Slope"),\
1701 (results['intercept'],results['eintercept'],\
1702 "Intercept")]:
1703 print("{0} = {1} +/- {2}".format(name, x, x_err))
1704
1705 corrcoef = np.corrcoef(allcoef[:,0], allcoef[:,1])
1706 results['corrcoef'] = corrcoef[1,0]
1707 results['covariance'] = results['corrcoef']*results['eslope']\
1708 *results['eintercept']
1709 results['results'] = allcoef
1710 print("The correlation coefficient between slope & intercept is {0}."\
1711 .format(results['corrcoef']))
1712 print("This leads to a covariance of {0} for slope & intercept."\
1713 .format(results['covariance']))
1714 print("Finally, {0} data points were available to produce this fit."\
1715 .format(results['n_points']))
1716 else:
1717 results['poly_degree'] = kwargs.get('poly_degree')
1718 results['results'] = allcoef
1719 return results
1720