;+ ;Language : IDL 8.2.2 ;================================================================================ ;The AAS gives permission to anyone who wishes to use these subroutines to run their own calculations. ;Permission to republish or reuse these routines should be directed to permissions@aas.org. ; ;Note that the AAS does not take responsibility for the content of the source code. Potential users should ;be wary of applying the code to conditions that the code was not written to model and the accuracy of the ;code may be affected when compiled and executed on different systems. ;================================================================================ ; NAME: ; BANYAN_MASS ; ; PURPOSE: ; This routine uses the J, H, Ks, W1 and W2 apparent magnitudes and the distance of a given object, along with ; an estimated age range to determine an 1-sigma mass range using AMES-COND isochrones (Baraffe et al. 2003) ; in combination with CIFIST2011 BT-SETTL atmosphere models (Allard et al. 2013, Rajpurohit et al. 2013) in a ; likelihood analysis. Errors on the distance and photometry, when input, are propagated to the estimated mass ; range. This routine was used to estimate the masses of candidate members in Table 3 of Gagné et al. (2014). ; ; CALLING SEQUENCE: ; BANYAN_MASS, age_0, age_1, dist, edist, mass_0, mass_1[, J=J, H=H, K=K, W1=W1, W2=W2, EJ=EJ, EH=EH, ; EK=EK, EW1=EW1, EW2=EW2, MAX_MASS=max_mass, BINARY=binary, S2M=s2m, SWISE=swise ; ; INPUTS: ; AGE_0 - The lower value of the estimated age range. ; AGE_1 - The higher value of the estimated age range. ; DIST - The distance measurement. ; EDIST - The error on the distance measurement. ; ; OPTIONAL INPUT KEYWORD: ; J - Apparent 2MASS J magnitude. ; H - Apparent 2MASS H magnitude. ; K - Apparent 2MASS Ks magnitude. ; W1 - Apparent WISE W1 magnitude. ; W2 - Apparent WISE W2 magnitude. ; EJ - Error on apparent 2MASS J magnitude. ; EH - Error on apparent 2MASS H magnitude. ; EK - Error on apparent 2MASS Ks magnitude. ; EW1 - Error on apparent WISE W1 magnitude. ; EW2 - Error on apparent WISE W2 magnitude. ; S2M - 2MASS structure containing photometry. ; SWISE - WISE structure containing photometry. ; MAX_MASS - The maximum mass that is tested on the grid, in Solar masses. By default this is set to 0.2 Solar masses. ; If a star with a larger mass is tested with this program, results will not be good. ; /BINARY - Set this keyword to remove 0.7 magnitudes to every channel if the tested object is an equal-luminosity binary. ; /BINARY - Set this keyword to receive output mass estimates in Solar masses instead of Jupiter masses. ; ; OUTPUTS: ; MASS_0 - The lower value of the estimated mass range, in Jupiter masses (unless /SOLAR is set). ; MASS_1 - The higher value of the estimated mass range, in Jupiter masses (unless /SOLAR is set). ; ; NOTE: ; By default, the Euler angles are expected to be input in radians. ; Please cite the following papers when using this procedure : ; Gagné et al. 2014, Baraffe et al. 2003, Allard et al. 2013, Rajpurohit et al. 2013 ; ; This routine requires the "barrafe_interpol_input.sav" IDL data save file. Please obtain this file at ; http://sites.google.com/site/mbderg/banyanii_additional_material under section 5.6 and update the ; "idl_datasave_path" variable at the start of "baraffe_interpol.pro" to a valid path on your computer ; leading to the file. ; ; PROCEDURES USED: ; DELETEVAR, BARAFFE_INTERPOL, INTEGRAL(), STDDEV() ; ; MODIFICATION HISTORY: ; WRITTEN, Jonathan Gagne, August, 21 2013 ;- Pro banyan_mass, age_0, age_1, dist, edist, mass_0, mass_1, J=J, H=H, K=K, W1=W1, W2=W2, EJ=EJ, EH=EH, $ EK=EK, EW1=EW1, EW2=EW2, MAX_MASS=max_mass, BINARY=binary, S2M=s2m, SWISE=swise ;Prend en entree des magnitudes absolues et des limites sur l'age, puis retourne un range de masses ;Un objet a la fois forward_function deletevar, baraffe_interpol, integral, stddev ;Define NaN value nan = !values.f_nan if ~keyword_set(dist) then $ message, 'Syntax : banyan_mass, age_0, age_1, dist, edist, mass_0, mass_1' if keyword_set(s2m) then begin J = s2m.J_M H = s2m.H_M K = s2m.K_M EJ = s2m.J_CMSIG EH = s2m.H_CMSIG EK = s2m.K_CMSIG endif if keyword_set(swise) then begin W1 = swise.W1MPRO W2 = swise.W2MPRO EW1 = swise.W1SIGMPRO EW2 = swise.W2SIGMPRO endif if ~keyword_set(j) and ~keyword_set(h) and ~keyword_set(k) and ~keyword_set(w1) and ~keyword_set(w2) then begin mass_0 = !values.f_nan mass_1 = !values.f_nan return endif dmod = 5.*alog10(dist)-5. edmod = edist/dist*5./alog(10) if J eq !NULL or EJ eq !NULL then begin J = nan EJ = nan endif else if keyword_set(dmod) then begin J -= dmod EJ = sqrt(EJ^2 + edmod^2) endif if H eq !NULL or EH eq !NULL then begin H = nan EH = nan endif else if keyword_set(dmod) then begin H -= dmod EH = sqrt(EH^2 + edmod^2) endif if K eq !NULL or EK eq !NULL then begin K = nan EK = nan endif else if keyword_set(dmod) then begin K -= dmod EK = sqrt(EK^2 + edmod^2) endif if W1 eq !NULL or EW1 eq !NULL then begin W1 = nan EW1 = nan endif else if keyword_set(dmod) then begin W1 -= dmod EW1 = sqrt(EW1^2 + edmod^2) endif if W2 eq !NULL or EW2 eq !NULL then begin W2 = nan EW2 = nan endif else if keyword_set(dmod) then begin W2 -= dmod EW2 = sqrt(EW2^2 + edmod^2) endif if keyword_set(binary) then begin J += alog10(2.)*2.5 H += alog10(2.)*2.5 K += alog10(2.)*2.5 W1 += alog10(2.)*2.5 W2 += alog10(2.)*2.5 endif ;Create an array of ages, for each of which the mass will be calculated np = 100. ages = findgen(np)/(float(np)-1.)*(max([age_0,age_1])-min([age_0,age_1]))+min([age_0,age_1]) masses = fltarr(np)+!values.f_nan emasses_pos = fltarr(np)+!values.f_nan emasses_neg = fltarr(np)+!values.f_nan mass_vec = !NULL ;Start loop to compute individual masses at each given age for i=0L, np-1L do begin deletevar, mass_vec, jt, ht, kt, w1t, w2t ;Interpolate the magnitudes VS mass relationships at given age, from AMES-COND isochrones (Baraffe et al. 2003) in combination with CIFIST2011 BT-SETTL atmosphere models (Allard et al. 2013, Rajpurohit et al. 2013) baraffe_interpol, ages[i], mass_vec, jt, ht, kt, w1t, w2t, Lit, MAX_MASS=max_mass if ~keyword_set(pdf_2d) then begin pdf_2d = dblarr(np,n_elements(mass_vec)) + !values.d_nan pdf_J = pdf_2d pdf_H = pdf_2d pdf_K = pdf_2d pdf_W1 = pdf_2d pdf_W2 = pdf_2d endif ;Build individual PDFs for each band pdf = dblarr(n_elements(jt))+0d0 OK = 0L if finite(J*EJ) then begin probJ = alog10(1./(sqrt(2.*!pi)*EJ)*exp(-(J-JT)^2/(2d0*EJ^2))) pdf += probJ OK = 1 endif if finite(H*EH) then begin probH = alog10(1./(sqrt(2.*!pi)*EH)*exp(-(H-HT)^2/(2d0*EH^2))) pdf += probH OK = 1 endif if finite(K*EK) then begin probK = alog10(1./(sqrt(2.*!pi)*EK)*exp(-(K-KT)^2/(2d0*EK^2))) pdf += probK OK = 1 endif if finite(W1*EW1) then begin probW1 = alog10(1./(sqrt(2.*!pi)*EW1)*exp(-(W1-W1T)^2/(2d0*EW1^2))) pdf += probW1 OK = 1 endif if finite(W2*EW2) then begin probW2 = alog10(1./(sqrt(2.*!pi)*EW2)*exp(-(W2-W2T)^2/(2d0*EW2^2))) pdf += probW2 OK = 1 endif ;Return NaN values if no magnitude can be used if OK eq 0L then begin mass_0 = !values.f_nan mass_1 = !values.f_nan return endif ;Compute cumulative density function and find most probable mass pdf = 10.^pdf g = where(finite(pdf)) cdf = integral(mass_vec[g], pdf[g], /accumulate) void = min(abs(cdf-max(cdf,/nan)/2d0),wm) masses[i] = mass_vec[g[wm]] ;Assign PDF for current age to a slice of the 2D PDF pdf_2d[i,*] = double(pdf) ;Do the same with individual bands if keyword_set(probJ) then $ pdf_J[i,*] = 10.^double(probJ) if keyword_set(probh) then $ pdf_H[i,*] = 10.^double(probH) if keyword_set(probk) then $ pdf_K[i,*] = 10.^double(probK) if keyword_set(probw1) then $ pdf_W1[i,*] = 10.^double(probW1) if keyword_set(probw2) then $ pdf_W2[i,*] = 10.^double(probW2) ;Find for which width 68% of the area under the PDF curve is covered (hence the 1-sigma error associated with the mass_vec estimation at this age) void = min((abs(abs(cdf-cdf[wm])/max(cdf,/nan)*2d0-erf(1d0/sqrt(2d0))))[0L:wm],wm1) void = min((abs(abs(cdf-cdf[wm])/max(cdf,/nan)*2d0-erf(1d0/sqrt(2d0))))[wm:*],wm2) wm2 += wm emasses_pos[i] = (mass_vec[g[wm]]-mass_vec[g[wm1]]) > (mass_vec[1]-mass_vec[0]) emasses_neg[i] = (mass_vec[g[wm2]]-mass_vec[g[wm]]) > (mass_vec[1]-mass_vec[0]) endfor ;Build prior PDFs for age age_prior = 1./(sqrt(2.*!pi)*abs(age_0-age_1))*exp(-(ages-mean([age_0,age_1]))^2/(2d0*abs(age_0-age_1)^2)) age_prior2d = age_prior#make_array((size(pdf_2d))[2],value=1.,/float) ;Marginalize PDF over age PDF_marg = total(pdf_2d*age_prior2d,1,/nan,/double) xv = mass_vec cdf = integral(xv, pdf_marg,/accumulate) nmes = max(cdf,/nan) ;Determine most probable mass where 50% of the area under curve of the PDF is covered void = min(abs(cdf-double(nmes)/2d0),wm) xf = xv[wm] ;Find for which width 68% of the area under the PDF curve is covered (hence the 1-sigma error associated with the final mass estimation) void = min((abs(abs(cdf-cdf[wm])/double(nmes)*2d0-erf(1d0/sqrt(2d0))))[0L:wm],wm1) void = min((abs(abs(cdf-cdf[wm])/double(nmes)*2d0-erf(1d0/sqrt(2d0))))[wm:*],wm2) wm2 += wm exf = [xv[wm]-xv[wm1],xv[wm2]-xv[wm]] mass_0 = xf - exf[0] mass_1 = xf + exf[1] if ~keyword_set(solar) then begin mass_0 *= 1047.2 mass_1 *= 1047.2 endif ;Avoid numerical errors when final masses are on the border of the grid if stddev(masses) lt median(masses)/1000. then begin mass_0 = !values.f_nan mass_1 = !values.f_nan return endif End