import math import os from astropy.coordinates import get_sun from astropy.time import Time from astropy.table import Table,Column, vstack import numpy as np import sys from time import time as xxtime import subprocess Folder = '/Users/jonas/' good_hpms_file = '/Users/jonas/hpms_good' RawCands_table = '/Users/jonas/rawpairs.vot' limit_shift = 0.03 #mas zeropoint = -0.029 DEG = math.pi/180e0 G_Newton = 6.67408e-11 # Gravitational konstant m^3 kg^-1 s^-2 M_sun = 1.98892e30 # Solar Mass in kg c_light = 299792458 # Speed of light in m/s pcinm = 3.08567758149137e16 # (1 parsec in m) m const_Einsteinradius = 180/math.pi * 3.6e6 * math.sqrt(4 * G_Newton * M_sun / (c_light ** 2 * pcinm * 1000)) #(mas/Msun)**0.5 def spherToCart(theta, phi): """returns a 3-cartesian unit vector pointing to longitude theta, latitude phi. The angles are in rad. """ cp = math.cos(phi) return math.cos(theta)*cp, math.sin(theta)*cp, math.sin(phi) def dirVecToCelCoos(dirVec): """returns alpha, delta in degrees for the direction vector dirVec. >>> dirVecToCelCoos(computeUnitSphereCoords(25.25, 12.125)) (25.25, 12.125) >>> dirVecToCelCoos(computeUnitSphereCoords(25.25, 12.125)*16) (25.25, 12.125) >>> "%g,%g"%dirVecToCelCoos(computeUnitSphereCoords(25.25, 12.125)+ ... computeUnitSphereCoords(30.75, 20.0)) '27.9455,16.0801' """ dirVec = dirVec.normalized() alpha = math.atan2(dirVec.y, dirVec.x) if alpha<0: alpha += 2*math.pi return alpha*1.8e2/math.pi, math.asin(dirVec.z)*1.8e2/math.pi class Vector3(object): """is a 3d vector that responds to both .x... and [0]... >>> x, y = Vector3(1,2,3), Vector3(2,3,4) >>> x+y Vector3(3.000000,5.000000,7.000000) >>> 4*x Vector3(4.000000,8.000000,12.000000) >>> x*4 Vector3(4.000000,8.000000,12.000000) >>> x*y 20 >>> "%.6f"%abs(x) '3.741657' >>> print abs((x+y).normalized()) 1.0 """ def __init__(self, x, y=None, z=None): if isinstance(x, tuple): self.coos = x else: self.coos = (x, y, z) def __repr__(self): return "Vector3(%f,%f,%f)"%tuple(self.coos) def __str__(self): def cutoff(c): if abs(c)<1e-10: return 0 else: return c rounded = [cutoff(c) for c in self.coos] return "[%.2g,%.2g,%.2g]"%tuple(rounded) def __getitem__(self, index): return self.coos[index] def __mul__(self, other): """does either scalar multiplication if other is not a Vector3, or a scalar product. """ if isinstance(other, Vector3): return self.x*other.x+self.y*other.y+self.z*other.z else: return Vector3(self.x*other, self.y*other, self.z*other) __rmul__ = __mul__ def __truediv__(self, scalar): return Vector3(self.x/scalar, self.y/scalar, self.z/scalar) def __add__(self, other): return Vector3(self.x+other.x, self.y+other.y, self.z+other.z) def __sub__(self, other): return Vector3(self.x-other.x, self.y-other.y, self.z-other.z) def __abs__(self): return math.sqrt(self.x**2+self.y**2+self.z**2) def cross(self, other): return Vector3(self.y*other.z-self.z*other.y, self.z*other.x-self.x*other.z, self.x*other.y-self.y*other.x) def normalized(self): return self/abs(self) def getx(self): return self.coos[0] def setx(self, x): self.coos[0] = x x = property(getx, setx) def gety(self): return self.coos[1] def sety(self, y): self.coos[1] = y y = property(gety, sety) def getz(self): return self.coos[2] def setz(self, z): self.coos[2] = z z = property(getz, setz) def calculate_Einsteinradius(lensMass, lensMass_err,lensparallax, lensparallax_err,obparallax = 0, obparallax_err = 0): """ calculate the Einstein radius in unit mas from Mass in SolarMasses, and paralax in mas ThetaE[rad] = sqrt(4GM/c^2)*sqrt(d_LS/(d_L*d_S)) ThetaE[mas] = const *sqrt(M/M_sun *(Lensparallax[mas]/1000-Objparallax[mas])) const = 180/pi * 3.6e6 *sqrt(4GM_sun / (c^2* 1pc[m]*1000)) """ if lensparallax < 0: ThetaE = const_Einsteinradius * (lensMass*(-lensparallax - obparallax)) ** 0.5 ThetaE_err = ThetaE * 0.5 * (lensMass_err ** 2 /lensMass ** 2 + (lensparallax_err ** 2 + obparallax_err ** 2) / (lensparallax - obparallax) **2 ) ** 0.5 else: ThetaE = const_Einsteinradius * (lensMass*abs(lensparallax - obparallax)) ** 0.5 ThetaE_err = ThetaE * 0.5 * (lensMass_err ** 2 /lensMass ** 2 + (lensparallax_err ** 2 + obparallax_err ** 2) / (lensparallax - obparallax) **2 ) ** 0.5 return ThetaE,ThetaE_err def calculate_Effekt(ThetaE, ThetaE_err, distance, distance_err,G,objG): """ calculate the microlensing effect for an given distance and Einsteinradius both in mas It return the distance in values of the Einstein radius, the shift in mas, the magnification in delta mag as well as their errors in the same units """ FL_FS = pow(100,(objG-G)/5) if distance is None: return (None,None,None,None,None,None,None,None,None,None) u = distance/ThetaE #distance in Einsteinradii du_u2 = (distance_err**2 / distance ** 2 + ThetaE_err ** 2 / ThetaE ** 2 ) # relativ error squard u_err = du_u2 ** 0.5 * u if u == 0: A = 1000 A_err = 1000 u_err = distance_err/ThetaE shift = math.sqrt(2)/4*ThetaE u_s = math.sqrt(2)/(1+FL_FS) # shift in mas shift_lum = u_s*ThetaE/(1+FL_FS)*(1+FL_FS*(u_s**2+3 - u_s * math.sqrt(u_s**2+4)))/(u_s**2+2+FL_FS*u_s*math.sqrt(u_s**2+4)) # shift in mas shift_p = ThetaE # shift shift_err= shift shift_lum_err = shift_lum shift_p_err=shift_p magnification = 1000 magnification_err = 1000 else: # absolut error if u < math.sqrt(2): # maximum shift at u = sqrt(2) shift = math.sqrt(2)/4*ThetaE # shift in mas shift_err = (((2 - u ** 2) / (u ** 2 + 2)) ** 2 * du_u2 + ThetaE_err ** 2 / ThetaE ** 2) ** 0.5 * shift else: shift = u / (u ** 2 + 2) * ThetaE # shift shift_err = (((2 - u ** 2) / (u ** 2 + 2)) ** 2 * du_u2 + ThetaE_err ** 2 / ThetaE ** 2) ** 0.5 * shift if u < math.sqrt(2)/(1+FL_FS): u_s = math.sqrt(2)/(1+FL_FS) shift_lum = u_s*ThetaE/(1+FL_FS)*(1+FL_FS*(u_s**2+3 - u_s * math.sqrt(u_s**2+4)))/(u_s**2+2+FL_FS*u_s*math.sqrt(u_s**2+4)) # shift in mas nn= ((ThetaE*FL_FS*u*(-u**2/math.sqrt(u**2 + 4) - math.sqrt(u**2 + 4) + 2*u))/((FL_FS + 1)*(FL_FS*math.sqrt(u**2 + 4)*u + u**2 + 2)) + (ThetaE*(FL_FS*(u**2 - math.sqrt(u**2 + 4)*u + 3) + 1))/((FL_FS + 1)*(FL_FS*math.sqrt(u**2 + 4)*u + u**2 + 2)) - (ThetaE*u*((FL_FS*u**2)/math.sqrt(u**2 + 4) + FL_FS*math.sqrt(u**2 + 4) + 2*u)*(FL_FS*(u**2 - math.sqrt(u**2 + 4)*u + 3) + 1)) /((FL_FS + 1)*(FL_FS*math.sqrt(u**2 + 4)*u + u**2 + 2)**2))**2 shift_lum_err = math.sqrt(nn*u_err**2+shift_lum**2*ThetaE_err**2/ThetaE**2) else: shift_lum = u*ThetaE/(1+FL_FS)*(1+FL_FS*(u**2+3 - u * math.sqrt(u**2+4)))/(u**2+2+FL_FS*u*math.sqrt(u**2+4)) # shift in mas nn= ((ThetaE*FL_FS*u*(-u**2/math.sqrt(u**2 + 4) - math.sqrt(u**2 + 4) + 2*u))/((FL_FS + 1)*(FL_FS*math.sqrt(u**2 + 4)*u + u**2 + 2)) + (ThetaE*(FL_FS*(u**2 - math.sqrt(u**2 + 4)*u + 3) + 1))/((FL_FS + 1)*(FL_FS*math.sqrt(u**2 + 4)*u + u**2 + 2)) - (ThetaE*u*((FL_FS*u**2)/math.sqrt(u**2 + 4) + FL_FS*math.sqrt(u**2 + 4) + 2*u)*(FL_FS*(u**2 - math.sqrt(u**2 + 4)*u + 3) + 1)) /((FL_FS + 1)*(FL_FS*math.sqrt(u**2 + 4)*u + u**2 + 2)**2))**2 shift_lum_err = math.sqrt(nn*u_err**2+shift_lum**2*ThetaE_err**2/ThetaE**2) shift_p = (math.sqrt(u**2+4)-u)/2*ThetaE shift_p_err = math.sqrt(u_err**2/(u**2+4)+ThetaE_err**2/ThetaE**2)* shift_p A = (u**2 +2)/(u*(u ** 2+4)**0.5) A_err = ((( -(u ** 2 + 2) / (u ** 2 * (u ** 2 + 4) ** 0.5) - (u ** 2 + 2) / (u ** 2 + 4) ** 1.5 + 2 / (u ** 2 + 4) ** 0.5) * u_err) ** 2) ** 0.5 magnification = 5 * math.log((FL_FS + A)/(FL_FS + 1)) / math.log(100) magnification_err = 5 / math.log(100) * A_err / (FL_FS+A) return u, u_err, shift, shift_err,shift_lum,shift_lum_err,shift_p, shift_p_err, magnification, magnification_err def f(x, a1,b1,c1,a2,b2,n): mass = pow(pow(np.exp(a1*x*x+b1*x +c1),-n)+pow(np.exp(a2*x+b2),-n), -1/n) mass = max(0.07, mass) return mass def calculate_Mass(g, bp,rp,parallax): """ estimate the mass of an lense from photometric properties classify star in White Dwarfs, Red Giant and Main sequence For MS determine approximatly Mass from magnitud. """ absmag = bp + 5 * math.log(parallax/100, 10) MG = g + 5 * math.log(parallax/100, 10) g_rp = g-rp if bp == 999: popt = [7.86232823e-03, -2.90891912e-01, 1.18248766e+00, -3.01175062e-01, 1.88952947e+00, 8.71127084e+01] Mass = f(MG, *popt) if Mass == 0.07: Mass_err = 0.03 Type = "BD" else: Mass_err = 0.1 * Mass Type = "MS" elif 4.*g_rp*g_rp+4.5*g_rp+ 6. < absmag : #white dwarfs if 17 < absmag : Mass = 0.07 Mass_err = 0.03 Type = "BD" else: Mass = 0.65 Mass_err = 0.15 Type = "WD" elif -3.*g_rp*g_rp+8.*g_rp - 1.3 > absmag: #red giant Mass = 1. Mass_err = 0.5 Type = "RG" else: #mainsequence popt = [7.86232823e-03, -2.90891912e-01, 1.18248766e+00, -3.01175062e-01, 1.88952947e+00, 8.71127084e+01] Mass = f(MG, *popt) if Mass == 0.07: Mass_err = 0.03 Type = "BD" else: Mass_err = 0.1 * Mass #10% Fehler? Type = "MS" return Mass, Mass_err, Type def cosdeg(arg): return math.cos(arg/180*math.pi) def sindeg(arg): return math.sin(arg/180*math.pi) def getGCDist_3Vec(pos1, pos2): """returns the distance along a great circle between two points. The distance is in degrees, the input positions are Vektor3. """ scalarprod = (pos1.normalized())*(pos2.normalized()) if scalarprod > 0.9: dif = pos1.normalized() - pos2.normalized() return 2e0*math.asin(abs(dif)/2e0)/DEG return math.acos(scalarprod)/DEG def movePm_3Vec(raDeg, decDeg, pmra, pmdec, timeDiff, foreshort = 0): """returns cartesian coordinates for an object with pos ra, dec and pm pmra after timeDiff. pmra has to have cos(dec) applied, position in deg prop.motion in mas , the time unit is yours to choose. """ ra, dec = raDeg/1.8e2*math.pi, decDeg/1.8e2*math.pi sd, cd = math.sin(dec), math.cos(dec) sa, ca = math.sin(ra), math.cos(ra) pmra, pmdec = pmra/1.8e2*math.pi/3.6e6 , pmdec/1.8e2*math.pi/3.6e6 muAbs = math.sqrt(pmra**2+pmdec**2); muTot = muAbs+0.5e0*foreshort*timeDiff; if muAbs<1e-20: dirVec = Vector3(cd*ca,cd*sa,sd) return dirVec # this is according to Mueller, 115 (4.94) dirA = pmra/muAbs; dirD = pmdec/muAbs; sinMot = math.sin(muTot*timeDiff); cosMot = math.cos(muTot*timeDiff); dirVec = Vector3(-sd*ca*dirD*sinMot - sa*dirA*sinMot + cd*ca*cosMot, -sd*sa*dirD*sinMot + ca*dirA*sinMot + cd*sa*cosMot, +cd*dirD*sinMot + sd*cosMot) return dirVec def parallax_correction(helioVec,earthVec,parallax): """ returns parrallax corrected cartesian Coordinates as Vector3 helioVec and earthVec as Vector3, parallax in mas """ baryVec = helioVec+earthVec*parallax/1e3 return baryVec def movePm_parallax(raDeg, decDeg, pmra, pmdec, parallax, time2015,foreshort = 0, gaia = False): helioVec = movePm_3Vec(raDeg, decDeg, pmra, pmdec, time2015, foreshort = foreshort) if parallax >= zeropoint: earthVec = pos_earth(time2015) if gaia == True: baryVec = parallax_correction(helioVec,earthVec,parallax*1.01) else: baryVec = parallax_correction(helioVec,earthVec,parallax) return baryVec else: return helioVec def pos_earth(time2015): #retruns earth positions for t in julian years (365.25 days) after J2015.5 #if time2015 > 85: print(time2015) #if time2015 < -15: print(time2015) earthcoord = get_sun(Time(time2015+2015.5, format = 'jyear')) return Vector3(spherToCart(earthcoord.ra.rad, earthcoord.dec.rad))*earthcoord.distance.pc def find_Minimum(f, left, right, minInterval = 3e-8): r = 0.61803399 c = 1-r x0 = left x3 = right x1 = c*x3 + r*x0 x2 = r*x3 + c*x0 f1 = f(x1) f2 = f(x2) n = 0 while math.fabs(x3-x0) > minInterval:#*(math.fabs(x1)+math.fabs(x2)): n += 1 if f2 < f1: x0 = x1 x1 = x2 x2 = c*x0 + r*x3 f1 = f2 f2 = f(x2) if f1 <= f2: x3 = x2 x2 = x1 x1 = c*x3 + r*x0 f2 = f1 f1 = f(x1) if f1 < f2: return x1 else: return x2 def estimate_errors_approx(row,minDate,minDist): if minDate is None: return None,None else: DeltaAlpha = float(row["ra"]-row["ob_ra"])*3.6e6 DeltaDelta = float(row["dec"]-row["ob_dec"])*3.6e6 cd, sd = cosdeg(float(row["dec"])), sindeg(float(row["dec"])) cd, sd = cosdeg(float(row["ra"])), sindeg(float(row["ra"])) if not row["ob_pmra"]: DeltaPmAlpha = float(row["pmra"]) DeltaPmDelta = float(row["pmdec"]) objpmalphaerr = 1e1 objpmdeltaerr = 1e1 else: DeltaPmAlpha = float(row["pmra"]-row["ob_pmra"]) DeltaPmDelta = float(row["pmdec"]-row["ob_pmdec"]) objpmalphaerr = float(row["ob_pmra_error"]) objpmdeltaerr = float(row["ob_pmdec_error"]) #mas/year lensalphaerr = float(row["ra_error"]) lensdeltaerr = float(row["dec_error"]) objalphaerr = float(row["ob_ra_error"]) objdeltaerr = float(row["ob_dec_error"]) lenspmalphaerr = float(row["pmra_error"]) lenspmdeltaerr = float(row["pmdec_error"]) if not row["ob_parallax"]: objparallax = 0.0e0 objparallaxerr = 0.0e0 else: objparallax = float(row["ob_parallax"]) objparallaxerr = float(row["ob_parallax_error"]) lensparallax = float(row["parallax"]) lensparallaxerr = float(row["parallax_error"]) DeltaPmTot = math.sqrt(pow(DeltaPmAlpha,2)+pow(DeltaPmDelta,2)) DeltaAlphaerr = math.sqrt(pow(lensalphaerr,2)+pow(objalphaerr,2) + pow((lensparallax-objparallax),2)) DeltaDeltaerr = math.sqrt(pow(lensdeltaerr,2)+pow(objdeltaerr,2) + pow((lensparallax-objparallax),2)) DeltaPmAlphaerr = math.sqrt(pow(lenspmalphaerr,2)+pow(objpmalphaerr,2)) DeltaPmDeltaerr = math.sqrt(pow(lenspmdeltaerr,2)+pow(objpmdeltaerr,2)) #minDate = (DeltaAlpha * DeltaPmAlpha + DeltaDelta * DeltaPmDelta) / DeltaPmTot ** 2 #minDist = (DeltaAlpha * DeltaPmDelta - DeltaDelta * DeltaPmAlpha) / DeltaPmTot minDateerr = math.sqrt( pow(pow(DeltaPmDelta,2) * DeltaAlpha - 2 * DeltaPmAlpha * DeltaPmDelta * DeltaDelta - pow(DeltaPmAlpha,2) * DeltaAlpha,2) / pow(DeltaPmTot,8) * pow(DeltaPmAlphaerr,2) + pow(pow(DeltaPmAlpha,2) * DeltaDelta - 2 * DeltaPmAlpha * DeltaPmDelta * DeltaAlpha - pow(DeltaPmDelta,2) * DeltaDelta,2) / pow(DeltaPmTot,8) * pow(DeltaPmDeltaerr,2) + pow(DeltaPmAlpha,2) / pow(DeltaPmTot,4) * pow(DeltaAlphaerr,2) + pow(DeltaPmDelta,2) / pow(DeltaPmTot,4) * pow(DeltaDeltaerr,2)) minDisterr = math.sqrt( pow(DeltaAlphaerr * DeltaPmDelta,2) / pow(DeltaPmTot,2) + pow(DeltaDeltaerr * DeltaPmAlpha,2) / pow(DeltaPmTot,2) + pow((DeltaDelta * pow(DeltaPmDelta,2) + DeltaAlpha * DeltaPmDelta * DeltaPmAlpha) * DeltaPmAlphaerr,2) / pow(DeltaPmTot,6) + pow((DeltaAlpha * pow(DeltaPmAlpha,2) + DeltaDelta * DeltaPmAlpha * DeltaPmDelta) * DeltaPmDeltaerr,2) / pow(DeltaPmTot,6)) return minDateerr, minDisterr def estimate_errors_parallax(row,minDate,minDist, delta_t_approx = 1/26.,gaia = False): if minDate is None: return None,None else: cd, sd = cosdeg(row["dec"]), sindeg(row["dec"]) ca, sa = cosdeg(row["ra"]), sindeg(row["ra"]) ocd, osd = cosdeg(row["ob_dec"]), sindeg(row["ob_dec"]) oca, osa = cosdeg(row["ob_ra"]), sindeg(row["ob_ra"]) objparallax = float(row["ob_parallax"] or zeropoint) lensparallax = float(row["parallax"]) objparallaxerr = float(row["ob_parallax_error"] or 2.) lensparallaxerr = float(row["parallax_error"]) objhelioVec = Vector3(ocd*oca,ocd*osa,osd) lenshelioVec = Vector3(cd*ca,cd*sa,sd) earthVec = pos_earth(minDate-2015.5) if gaia: lensra,lensdec = dirVecToCelCoos(parallax_correction(lenshelioVec,earthVec,lensparallax*1.01)) objra,objdec = dirVecToCelCoos(parallax_correction(objhelioVec,earthVec,objparallax*1.01)) else: lensra,lensdec = dirVecToCelCoos(parallax_correction(lenshelioVec,earthVec,lensparallax)) objra,objdec = dirVecToCelCoos(parallax_correction(objhelioVec,earthVec,objparallax)) DeltaAlpha = float(lensra-objra)*3.6e6*cd DeltaDelta = float(lensdec-objdec)*3.6e6 if row["ob_pmra"]: DeltaPmAlpha = float(row["pmra"]-row["ob_pmra"]) DeltaPmDelta = float(row["pmdec"]-row["ob_pmdec"]) lensalphaerr = float(row["ra_error"]) lensdeltaerr = float(row["dec_error"]) objalphaerr = float(row["ob_ra_error"]) objdeltaerr = float(row["ob_dec_error"]) lenspmalphaerr = float(row["pmra_error"]) lenspmdeltaerr = float(row["pmdec_error"]) objpmalphaerr = float(row["ob_pmra_error"]) objpmdeltaerr = float(row["ob_pmdec_error"]) DeltaPmTot = math.sqrt(pow(DeltaPmAlpha,2)+pow(DeltaPmDelta,2)) DeltaAlphaerr = math.sqrt(pow(lensalphaerr,2)+pow(objalphaerr,2)) DeltaDeltaerr = math.sqrt(pow(lensdeltaerr,2)+pow(objdeltaerr,2)) DeltaPmAlphaerr = math.sqrt(pow(lenspmalphaerr,2)+pow(objpmalphaerr,2)) DeltaPmDeltaerr = math.sqrt(pow(lenspmdeltaerr,2)+pow(objpmdeltaerr,2)) else: DeltaPmAlpha = float(row["pmra"]) DeltaPmDelta = float(row["pmdec"]) lensalphaerr = float(row["ra_error"]) lensdeltaerr = float(row["dec_error"]) objalphaerr = float(row["ob_ra_error"]) objdeltaerr = float(row["ob_dec_error"]) lenspmalphaerr = float(row["pmra_error"]) lenspmdeltaerr = float(row["pmdec_error"]) objpmalphaerr = 10e0 objpmdeltaerr = 10e0 DeltaPmTot = math.sqrt(pow(DeltaPmAlpha,2)+pow(DeltaPmDelta,2)) DeltaAlphaerr = math.sqrt(pow(lensalphaerr,2)+pow(objalphaerr,2)) DeltaDeltaerr = math.sqrt(pow(lensdeltaerr,2)+pow(objdeltaerr,2)) DeltaPmAlphaerr = math.sqrt(pow(lenspmalphaerr,2)+pow(objpmalphaerr,2)) DeltaPmDeltaerr = math.sqrt(pow(lenspmdeltaerr,2)+pow(objpmdeltaerr,2)) minDateerr = math.sqrt( pow(pow(DeltaPmDelta,2) * DeltaAlpha - 2 * DeltaPmAlpha * DeltaPmDelta * DeltaDelta - pow(DeltaPmAlpha,2) * DeltaAlpha,2) / pow(DeltaPmTot,8) * pow(DeltaPmAlphaerr,2) + pow(pow(DeltaPmAlpha,2) * DeltaDelta - 2 * DeltaPmAlpha * DeltaPmDelta * DeltaAlpha - pow(DeltaPmDelta,2) * DeltaDelta,2) / pow(DeltaPmTot,8) * pow(DeltaPmDeltaerr,2) + pow(DeltaPmAlpha,2) / pow(DeltaPmTot,4) * pow(DeltaAlphaerr,2) + pow(DeltaPmDelta,2) / pow(DeltaPmTot,4) * pow(DeltaDeltaerr,2)) dra,ddec = (lensra-row["ra"])*3.6e6*cd,(lensdec -row["dec"])*3.6e6 if minDateerr > 0.25: earthVec2 = pos_earth(minDate-2015.5+0.25) else: earthVec2 = pos_earth(minDate-2015.5+minDateerr) if gaia: lensra2,lensdec2 = dirVecToCelCoos(parallax_correction(lenshelioVec,earthVec2,(lensparallax-objparallax)*1.01)) #deg else: lensra2,lensdec2 = dirVecToCelCoos(parallax_correction(lenshelioVec,earthVec2,(lensparallax-objparallax))) dvra,dvdec = (lensra-lensra2)*3.6e6*cd ,(lensdec -lensdec2)*3.6e6 #deg parallaxerr2 = lensparallaxerr**2/lensparallax**2 + objparallaxerr**2/lensparallax**2 DeltaAlphaerr = math.sqrt(DeltaAlphaerr**2+ dra**2 * parallaxerr2 + dvra ** 2) DeltaDeltaerr = math.sqrt(DeltaDeltaerr**2 + ddec**2 * parallaxerr2 + dvdec ** 2) #minDate = (DeltaAlpha * DeltaPmAlpha + DeltaDelta * DeltaPmDelta) / DeltaPmTot ** 2 #minDist = (DeltaAlpha * DeltaPmDelta - DeltaDelta * DeltaPmAlpha ) / DeltaPmTot minDateerr = math.sqrt( pow(pow(DeltaPmDelta,2) * DeltaAlpha - 2 * DeltaPmAlpha * DeltaPmDelta * DeltaDelta - pow(DeltaPmAlpha,2) * DeltaAlpha,2) / pow(DeltaPmTot,8) * pow(DeltaPmAlphaerr,2) + pow(pow(DeltaPmAlpha,2) * DeltaDelta - 2 * DeltaPmAlpha * DeltaPmDelta * DeltaAlpha - pow(DeltaPmDelta,2) * DeltaDelta,2) / pow(DeltaPmTot,8) * pow(DeltaPmDeltaerr,2) + pow(DeltaPmAlpha,2) / pow(DeltaPmTot,4) * pow(DeltaAlphaerr,2) + pow(DeltaPmDelta,2) / pow(DeltaPmTot,4) * pow(DeltaDeltaerr,2)) minDisterr = math.sqrt( pow(DeltaAlphaerr * DeltaPmDelta,2) / pow(DeltaPmTot,2) + pow(DeltaDeltaerr * DeltaPmAlpha,2) / pow(DeltaPmTot,2) + pow((DeltaDelta * pow(DeltaPmDelta,2) + DeltaAlpha * DeltaPmDelta * DeltaPmAlpha) * DeltaPmAlphaerr,2) / pow(DeltaPmTot,6) + pow((DeltaAlpha * pow(DeltaPmAlpha,2) + DeltaDelta * DeltaPmAlpha * DeltaPmDelta) * DeltaPmDeltaerr,2) / pow(DeltaPmTot,6)) return minDateerr, minDisterr def estimate_Closest_approx(lensra, lensdec, lenspmra, lenspmdec, objra, objdec, objpmra, objpmdec): dist1 = lambda t: getGCDist_3Vec( movePm_3Vec(lensra, lensdec, lenspmra, lenspmdec, t), movePm_3Vec(objra, objdec, objpmra, objpmdec, t)) * 3.6e6 t_0 = find_Minimum(dist1, -15.002, 85.002, minInterval = 0.0001) approx_closest_date = t_0 + 2015.5 approx_closest_dist = dist1(t_0) return approx_closest_date, approx_closest_dist def estimate_Closest_parallax(lensra, lensdec, lenspmra, lenspmdec, lensparallax, objra, objdec, objpmra, objpmdec, objparallax,t_0,gaia = True): """ units: Ra,DEC in Deg ra_err, Dec_err in mas pmRA, pmDec in mas/year, pmRa_err,pmDec_err in mas/year parrallax in mas parrallax_err in mas distance in mas date in year mass in solar_mass einsteinradii in mas shift in mas magnifikation in delta mag """ # Define distance function with and with out parallax dist2 = lambda t: getGCDist_3Vec( movePm_parallax(lensra, lensdec, lenspmra, lenspmdec, lensparallax, t), movePm_parallax(objra, objdec, objpmra, objpmdec, objparallax, t)) *3.6e6 Mu_tot = ((lenspmra-objpmra)**2 * cosdeg(lensdec)**2 + (lenspmdec - objpmdec)**2)**0.5 #finde minima limit to approx 2 week if t_0 <= 85 and t_0 >= -15: if lensparallax != 0: tlim = max(1.1,(2*lensparallax) / Mu_tot) t_min = t_0 - tlim t_max = t_0 + tlim ti = list(map(lambda x: float(x/26.) ,range(int(t_min*26),int(t_max*26)))) dist = list(map(dist2,ti)) number_list = range(len(ti)-2) a = list(filter(lambda x: dist[x+1] - dist[x] < 0 and dist[x+1] - dist[x+2] <0 , number_list)) closest_time_vec = [999999.,999999.,999999.] closest_dist_vec = [999999.,999999.,999999.] for j in a: t1 = ti[j+0] t2 = ti[j+2] closest_time = find_Minimum(dist2,t1,t2) closest_dist = dist2(closest_time) if closest_dist <= closest_dist_vec[0]: closest_time_vec[2] = closest_time_vec[1] closest_time_vec[1] = closest_time_vec[0] closest_time_vec[0] = closest_time closest_dist_vec[2] = closest_dist_vec[1] closest_dist_vec[1] = closest_dist_vec[0] closest_dist_vec[0] = closest_dist elif closest_dist <= closest_dist_vec[1]: closest_time_vec[2] = closest_time_vec[1] closest_time_vec[1] = closest_time closest_dist_vec[2] = closest_dist_vec[1] closest_dist_vec[1] = closest_dist elif closest_dist <= closest_dist_vec[2]: closest_time_vec[2] = closest_time closest_dist_vec[2] = closest_dist closest_dist_vec = [x if x < 999999 else None for x in closest_dist_vec] closest_date_vec = [x+2015.5 if x < 999999 else None for x in closest_time_vec] if gaia == True: """ calculate closest approxhes seen from L2 therefor the parallax is multiplied with 1.01 """ dist3 = lambda t: getGCDist_3Vec( movePm_parallax(lensra, lensdec, lenspmra, lenspmdec, lensparallax, t, gaia = True), movePm_parallax(objra, objdec, objpmra, objpmdec, objparallax, t, gaia = True)) * 3.6e6 dist = list(map(dist3,ti)) number_list = range(len(ti)-2) a = list(filter(lambda x: dist[x+1] - dist[x] < 0 and dist[x+1] - dist[x+2] <0 , number_list)) closest_time_gaia_vec = [999999.,999999.,999999.] closest_dist_gaia_vec = [999999.,999999.,999999.] for j in a: t1 = ti[j+0] t2 = ti[j+2] closest_time = find_Minimum(dist3,t1,t2) closest_dist = dist3(closest_time) if closest_dist <= closest_dist_gaia_vec[0]: closest_time_gaia_vec[2] = closest_time_gaia_vec[1] closest_time_gaia_vec[1] = closest_time_gaia_vec[0] closest_time_gaia_vec[0] = closest_time closest_dist_gaia_vec[2] = closest_dist_gaia_vec[1] closest_dist_gaia_vec[1] = closest_dist_gaia_vec[0] closest_dist_gaia_vec[0] = closest_dist elif closest_dist <= closest_dist_gaia_vec[1]: closest_time_gaia_vec[2] = closest_time_gaia_vec[1] closest_time_gaia_vec[1] = closest_time closest_dist_gaia_vec[2] = closest_dist_gaia_vec[1] closest_dist_gaia_vec[1] = closest_dist elif closest_dist <= closest_dist_gaia_vec[2]: closest_time_gaia_vec[2] = closest_time closest_dist_gaia_vec[2] = closest_dist closest_dist_gaia_vec = [x if x < 999999 else None for x in closest_dist_gaia_vec] closest_date_gaia_vec = [x+2015.5 if x < 999999 else None for x in closest_time_gaia_vec] return closest_date_vec, closest_dist_vec, closest_date_gaia_vec, closest_dist_gaia_vec else: return closest_date_vec,closest_dist_vec else: if gaia: closest_time_vec = [None,None,None] closest_dist_vec = [None,None,None] closest_time_gaia_vec = [None,None,None] closest_dist_gaia_vec = [None,None,None] return closest_date_vec, closest_dist_vec, closest_date_gaia_vec, closest_dist_gaia_vec else: closest_time_vec = [None,None,None] closest_dist_vec = [None,None,None] return closest_date_vec,closest_dist_vec else: return () def find_Closest(row, lensra, lensdec, objra, objdec, lenspmra, lenspmdec, objpmra, objpmdec, lensparallax, objparallax,Gmag,objGmag ,gaia = True, ThetaE = None,ThetaEerr = None): """ units: Ra,DEC in Deg ra_err, Dec_err in mas pmRA, pmDec in mas/year, pmRa_err,pmDec_err in mas/year parrallax in mas parrallax_err in mas distance in mas date in year mass in solar_mass einsteinradii in mas shift in mas magnifikation in delta mag """ Taa, Daa = estimate_Closest_approx(lensra, lensdec, lenspmra, lenspmdec, objra, objdec, objpmra, objpmdec) Taaerr, Daaerr = estimate_errors_approx(row, Daa, Taa) T_0 = Taa - 2015 approxdE, approxdEerr, approxshift, approxshifterr,approxshiftlum, approxshiftlumerr,approxshiftp, approxshiftperr, approxmagn, approxmagnerr = calculate_Effekt(ThetaE, ThetaEerr, Daa, Daaerr,Gmag,objGmag) if approxdE < (ThetaE/0.1): timescale = 2 * np.sqrt((pow(ThetaE/0.1,2)-pow(approxdE,2))/(pow(lenspmra,2)+pow(lenspmdec,2)))*ThetaE else: timescale = 0 if T_0 <= -15 or T_0 > 70: if gaia: return (None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None) else: return (None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None) if (lensparallax*1.01 > approxdE or approxshift * approxdE/max(approxdE-lensparallax*1.01,1) > limit_shift ) and lensparallax > 0: cc = estimate_Closest_parallax( lensra, lensdec, lenspmra, lenspmdec, lensparallax, objra, objdec, objpmra, objpmdec, objparallax, T_0, gaia) if len(cc) == 0: if gaia: return (None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None) else: return (None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None) if gaia: minDatevec, minDistvec, minDategaiavec, minDistgaiavec = cc Tca, Tfa1, Tfa2 = minDatevec Dca, Dfa1, Dfa2 = minDistvec TcaGaia, Tfa1gaia, Tfa2Gaia = minDategaiavec DcaGaia, Dfa1Gaia, Dfa2Gaia = minDistgaiavec Tcaerr, Dcaerr = estimate_errors_parallax(row,Tca,Dca, delta_t_approx = Taaerr) TcaGaiaerr, DcaGaiaerr = estimate_errors_parallax(row,TcaGaia,DcaGaia, delta_t_approx = Taaerr, gaia = True) if DcaGaiaerr is None and DcaGaia is not None: print(Taa,Daa) print(lensra,lensdec,lenspmra, lenspmdec,lensparallax) print( objra, objdec,objpmra, objpmdec, objparallax) print(minDategaiavec) print(minDistgaiavec) #Dcaerr, Tcaerr, DcaGaiaerr, TcaGaiaerr = estimate_errors_parallax(row, Taa, gaia) DcaE, DcaEerr, Sca, Scaerr,ScaLum,ScaLumerr,ScaP,ScaPerr, Mca, Mcaerr = calculate_Effekt(ThetaE, ThetaEerr, Dca, Dcaerr,Gmag,objGmag) DcaEGaia, DcaEGaiaerr, ScaGaia, ScaGaiaerr,ScaLumGaia, ScaLumGaiaerr,ScaPGaia, ScaPGaiaerr, McaGaia, McaGaiaerr = calculate_Effekt(ThetaE, ThetaEerr, DcaGaia, DcaGaiaerr,Gmag,objGmag) Dfa1E, Dfa1Eerr, Sfa1, Sfa1err,Sfa1Lum, Sfa1Lumerr,Sfa1P, Sfa1Perr, Mfa1, Mfa1err = calculate_Effekt(ThetaE, ThetaEerr, Dfa1,Daaerr,Gmag,objGmag) Dfa1EGaia, Dfa1EGaiaerr, Sfa1Gaia, Sfa1Gaiaerr,Sfa1LumGaia, Sfa1LumGaiaerr,Sfa1PGaia, Sfa1PGaiaerr, Mfa1Gaia, Mfa1Gaiaerr = calculate_Effekt(ThetaE, ThetaEerr, Dfa1Gaia, DcaGaiaerr,Gmag,objGmag) if Sfa1 is None: Fa = False elif Sfa1 >= limit_shift: Fa = True ''' write row to file ''' else: Fa = False if Sfa1Gaia is None: FaGaia = False elif Sfa1Gaia >= limit_shift : FaGaia = True ''' write row to file ''' else: FaGaia = False return (Taa,Daa,Taaerr,Daaerr, Tca, Dca, Tcaerr, Dcaerr, Fa, TcaGaia, DcaGaia, TcaGaiaerr, DcaGaiaerr,FaGaia, approxdE, approxdEerr, approxshift, approxshifterr,approxshiftlum, approxshiftlumerr,approxshiftp, approxshiftperr, approxmagn, approxmagnerr, DcaE, DcaEerr, Sca, Scaerr,ScaLum, ScaLumerr,ScaP, ScaPerr, Mca, Mcaerr, DcaEGaia, DcaEGaiaerr, ScaGaia, ScaGaiaerr, ScaLumGaia, ScaLumGaiaerr,ScaPGaia, ScaPGaiaerr, McaGaia, McaGaiaerr, timescale) else: minDatevec, minDistvec = cc Tca, Tfa1, Tfa2 = minDatevec Dca, Dfa1, Dfa2 = minDistvec Tcaerr, Dcaerr = estimate_errors_parallax(row,Tca,Dca, delta_t_approx = Taaerr) DcaE, DcaEerr, Sca, Scaerr,ScaLum, ScaLumerr,ScaP, ScaPerr, Mca, Mcaerr = calculate_Effekt(ThetaE, ThetaEerr, Dca, Dcaerr,Gmag,objGmag) Dfa1E, Dfa1Eerr, Sfa1, Sfa1err,Sfa1Lum, Sfa1Lumerr,Sfa1P, Sfa1Perr, Mfa1, Mfa1err = calculate_Effekt(ThetaE, ThetaEerr, Dfa1,Daaerr,Gmag,objGmag) if Sfa1 is None: Sfa1 = 0.5 * limit_shift if Sfa1 >= 3*limit_shift: Fa = True ''' write row to file ''' else: Fa = False return (Taa,Daa,Taaerr,Daaerr, Tca, Dca, Tcaerr, Dcaerr, Fa, approxdE, approxdEerr, approxshift, approxshifterr,approxshiftlum, approxshiftlumerr,approxshiftp, approxshiftperr, approxmagn, approxmagnerr, DcaE, DcaEerr, Sca, Scaerr,ScaLum, ScaLumerr,ScaP, ScaPerr, Mca, Mcaerr) else: if gaia: Tca, Dca, Tcaerr, Dcaerr = [None,None,None,None] TcaGaia, DcaGaia, TcaGaiaerr, DcaGaiaerr = [None,None,None,None] DcaE, DcaEerr, Sca, Scaerr,ScaLum, ScaLumerr,ScaP, ScaPerr, Mca, Mcaerr = [None,None,None,None,None,None,None,None,None,None] DcaEGaia, DcaEGaiaerr, ScaGaia, ScaGaiaerr,ScaLumGaia, ScaLumGaiaerr,ScaPGaia, ScaPGaiaerr, McaGaia, McaGaiaerr = [None,None,None,None,None,None,None,None,None,None] FaGaia = False Fa = False return (Taa,Daa,Taaerr,Daaerr, Tca, Dca, Tcaerr, Dcaerr, Fa, TcaGaia, DcaGaia, TcaGaiaerr, DcaGaiaerr,FaGaia, approxdE, approxdEerr, approxshift, approxshifterr,approxshiftlum, approxshiftlumerr,approxshiftp, approxshiftperr, approxmagn, approxmagnerr, DcaE, DcaEerr, Sca, Scaerr, ScaLum, ScaLumerr,ScaP, ScaPerr, Mca, Mcaerr, DcaEGaia, DcaEGaiaerr, ScaGaia, ScaGaiaerr,ScaLumGaia, ScaLumGaiaerr,ScaPGaia, ScaPGaiaerr, McaGaia, McaGaiaerr, timescale) else: Tca, Dca, Tcaerr, Dcaerr = [None,None,None,None] DcaE, DcaEerr, Sca, Scaerr, ScaLum, ScaLumerr,ScaP, ScaPerr, Mca, Mcaerr = [None,None,None,None,None,None,None,None,None,None] Fa = False return (Taa,Daa,Taaerr,Daaerr, Tca, Dca, Tcaerr, Dcaerr, Fa, approxdE, approxdEerr, approxshift, approxshifterr,approxshiftlum, approxshiftlumerr,approxshiftp, approxshiftperr, approxmagn, approxmagnerr, DcaE, DcaEerr, Sca, Scaerr, ScaLum, ScaLumerr,ScaP, ScaPerr, Mca, Mcaerr,timescale) def iterrow(row): srcId = row["source_id"] lensra, lensdec = float(row["ra"]), float(row["dec"]) #deg lensraerr, lensdecerr = float(row["ra_error"]), float(row["dec_error"]) #mas lenspmra, lenspmdec = float(row["pmra"]), float(row["pmdec"]) #mas/year lenspmraerr, lenspmdecerr = float(row["pmra_error"]), float(row["pmdec_error"]) #mas/year lensparallax, lensparallaxerr = float(row["parallax"]), float(row["parallax_error"]) #mas lensG= float(row["phot_g_mean_mag"]) #mag lensBp, lensRp = float(row["phot_bp_mean_mag"] or 999.), float(row["phot_rp_mean_mag"] or 999.) objsrcId = row["ob_source_id"] objra, objdec = float(row["ob_ra"]), float(row["ob_dec"]) #deg objraerr, objdecerr = float(row["ob_ra_error"]), float(row["ob_dec_error"]) #mas objpmra, objpmdec = float(row["ob_pmra"] or 0.0), float(row["ob_pmdec"] or 0.0) objpmraerr, objpmdecerr = float(row["ob_pmra_error"] or 10.0), float(row["ob_pmdec_error"] or 10.0) #mas/year objparallax, objparallaxerr = float(row["ob_parallax"] or zeropoint), float(row["ob_parallax_error"] or 2.) #mas if objparallax < zeropoint: objparallax, zeropoint #mas objG = float(row["ob_phot_g_mean_mag"]) #mag objBp, objRp = float(row["ob_phot_bp_mean_mag"] or 999.), float(row["ob_phot_rp_mean_mag"] or 999.) #mag if lensparallax < objparallax: return None else: lensMass, lensMasserr, startype = calculate_Mass(lensG, lensBp, lensRp,lensparallax) if objparallax <= zeropoint: ThetaE,ThetaEerr = calculate_Einsteinradius(lensMass, lensMasserr,lensparallax+ zeropoint, lensparallaxerr) else: ThetaE,ThetaEerr = calculate_Einsteinradius(lensMass, lensMasserr,lensparallax, lensparallaxerr,objparallax, objparallaxerr) [approxDate,approxDist,approxDateerr,approxDisterr, Tca, Dca, Tcaerr, Dcaerr, Fa, TcaGaia, DcaGaia, TcaGaiaerr, DcaGaiaerr,FaGaia, approxdE, approxdEerr, approxshift, approxshifterr,approxshiftlum, approxshiftlumerr,approxshiftp, approxshiftperr, approxmagn, approxmagnerr, DcaE, DcaEerr, Sca, Scaerr, ScaLum, ScaLumerr,ScaP, ScaPerr, Mca, Mcaerr, DcaEGaia, DcaEGaiaerr, ScaGaia, ScaGaiaerr,ScaLumGaia, ScaLumGaiaerr,ScaPGaia, ScaPGaiaerr, McaGaia, McaGaiaerr, timescale ] = find_Closest(row, lensra, lensdec, objra, objdec, lenspmra, lenspmdec, objpmra, objpmdec, lensparallax, objparallax,lensG,objG, gaia = True, ThetaE = ThetaE, ThetaEerr = ThetaEerr) if approxDate is None: return None if approxDate < 2000.5 or approxDate > 2070: return None if ScaP is None: return None if ScaP < 0.1: return None row_out = Table(row) row_out['star_type'] = Column([startype], dtype = 'str', description = 'Type of the lensing star: WD = White Dwarf, MS = Main Sequence, RG = Red Giant, BD = Brown Dwarf') row_out['mass'] = Column([lensMass],dtype = 'float64', unit = 'M_sun', description = 'Mass of the lensing star') row_out['mass_error'] = Column([lensMasserr],dtype = 'float64', unit = 'M_sun', description = 'Error of the Mass of the lensing star') row_out['thetaE'] = Column([ThetaE],dtype = 'float64', unit = 'mas', description = 'Einstein radius of the event') row_out['thetaE_error'] = Column([ThetaEerr],dtype = 'float64', unit = 'mas', description = 'Error of the Einstein radius') row_out['time_scale'] = Column([timescale], dtype = 'float64' , unit = 'year', description = 'Approximated duration of an event') row_out['approx_date'] = Column([approxDate], dtype = 'float64' , unit = 'year', description = 'Time of the approximated closest approch') row_out['approx_date_error'] = Column([approxDateerr],dtype = 'float64', unit = 'year', description = 'Error of the Time of the approximated closest approch') row_out['approx_dist'] = Column([approxDist],dtype = 'float64', unit = 'mas', description = 'Approximated closest distance') row_out['approx_dist_error'] = Column([approxDisterr],dtype = 'float64', unit = 'mas', description = 'Error of approximated closest distance') row_out['approx_u'] = Column([approxdE],dtype = 'float64', unit = '', description = 'Approximated closest distance in einstein radii') row_out['approx_u_error'] = Column([approxdEerr],dtype = 'float64', unit = '', description = 'Error of the approximated closest distance in einstein radii') row_out['approx_shift'] = Column([approxshift],dtype = 'float64', unit = 'mas', description = 'Approximated maximal astrometric shift of the center of light') row_out['approx_shift_error'] = Column([approxshifterr],dtype = 'float64', unit = 'mas', description = 'Error of the approximated shift of the center of light') row_out['approx_shift_lum'] = Column([approxshiftlum],dtype = 'float64', unit = 'mas', description = 'Approximated maximal astrometric shift includig lens-luminosity effects') row_out['approx_shift_lum_error'] = Column([approxshiftlumerr],dtype = 'float64', unit = 'mas', description = 'Error of the approximated shift includig lens-luminosity effects') row_out['approx_shift_plus'] = Column([approxshiftp],dtype = 'float64', unit = 'mas', description = 'Approximated maximal astrometric shift of image (+)') row_out['approx_shift_plus_error'] = Column([approxshiftperr],dtype = 'float64', unit = 'mas', description = 'Error of the approximated shift of image (+)') row_out['approx_magnification'] = Column([approxmagn],dtype = 'float64', unit = 'mag', description = 'Approximated magnification in magnitude') row_out['approx_magnification_error'] = Column([approxmagnerr],dtype = 'float64', unit = 'mag', description = 'Error of the approximated magnification') row_out['date'] = Column([Tca],dtype = 'float64', unit = 'year', description = 'Time of the closest approch') row_out['date_error'] = Column([Tcaerr],dtype = 'float64', unit = 'year', description = 'Error of the Time of theclosest approch') row_out['dist'] = Column([Dca],dtype = 'float64', unit = 'mas', description = 'Closest distance') row_out['dist_error'] = Column([Dcaerr],dtype = 'float64', unit = 'mas', description = 'Error of the closest distance') row_out['u'] = Column([DcaE],dtype = 'float64', unit = '', description = 'Closest distance in einstein radii') row_out['u_error'] = Column([DcaEerr],dtype = 'float64', unit = '', description = 'Error of the closest distance in einstein radii') row_out['shift'] = Column([Sca],dtype = 'float64', unit = 'mas', description = 'Expected maximal astrometric shift of the center of light') row_out['shift_error'] = Column([Scaerr],dtype = 'float64', unit = 'mas', description = 'Error of the shift of the center of light') row_out['shift_lum'] = Column([ScaLum],dtype = 'float64', unit = 'mas', description = 'Expected maximal astrometric shift includig lens-luminosity effects') row_out['shift_lum_error'] = Column([ScaLumerr],dtype = 'float64', unit = 'mas', description = 'Error of the shift includig lens-luminosity effects') row_out['shift_plus'] = Column([ScaP],dtype = 'float64', unit = 'mas', description = 'Expected maximal astrometric shift of image (+)') row_out['shift_plus_error'] = Column([ScaPerr],dtype = 'float64', unit = 'mas', description = 'Error of the shift of image (+)') row_out['magnification'] = Column([Mca],dtype = 'float64', unit = 'mag', description = 'Magnification in magnitudes') row_out['magnification_error'] = Column([Mcaerr],dtype = 'float64', unit = 'mag', description = 'Error of the magnification in magnitudes') row_out['L2_date'] = Column([TcaGaia],dtype = 'float64', unit = 'year', description = 'Time of the closest approch for Lagrange Point 2 parallax') row_out['L2_date_error'] = Column([TcaGaiaerr],dtype = 'float64', unit = 'year', description = 'Error of the Time of theclosest approch for Lagrange Point 2 parallax') row_out['L2_dist'] = Column([DcaGaia],dtype = 'float64', unit = 'mas', description = 'Closest distance for Lagrange Point 2 parallax') row_out['L2_dist_error'] = Column([DcaGaiaerr],dtype = 'float64', unit = 'mas', description = 'Error of the closest distance for Lagrange Point 2 parallax') row_out['L2_u'] = Column([DcaEGaia],dtype = 'float64', unit = '', description = 'Closest distance in einstein radii for Lagrange Point 2 parallax') row_out['L2_u_error'] = Column([DcaEGaiaerr],dtype = 'float64', unit = '', description = 'Error of the closest distance in einstein radii for Lagrange Point 2 parallax') row_out['L2_shift'] = Column([ScaGaia],dtype = 'float64', unit = 'mas', description = 'Expected maximal astrometric shift of the center of light for Lagrange Point 2 parallax') row_out['L2_shift_error'] = Column([ScaGaiaerr],dtype = 'float64', unit = 'mas', description = 'Error of the shift of the center of light for Lagrange Point 2 parallax') row_out['L2_shift_lum'] = Column([ScaLumGaia],dtype = 'float64', unit = 'mas', description = 'Expected maximal shift includig lens-luminosity effects for Lagrange Point 2 parallax') row_out['L2_shift_lum_error'] = Column([ScaLumGaiaerr],dtype = 'float64', unit = 'mas', description = 'Error of the shift includig lens-luminosity effects for Lagrange Point 2 parallax') row_out['L2_shift_plus'] = Column([ScaPGaia],dtype = 'float64', unit = 'mas', description = 'Expected maximal shift of image (+) for Lagrange Point 2 parallax') row_out['L2_shift_plus_error'] = Column([ScaPGaiaerr],dtype = 'float64', unit = 'mas', description = 'Error of the shift of image (+) for Lagrange Point 2 parallax') row_out['L2_magnification'] = Column([McaGaia],dtype = 'float64', unit = 'mag', description = 'Magnification in magnitudes for Lagrange Point 2 parallax') row_out['L2_magnification_error'] = Column([McaGaiaerr],dtype = 'float64', unit = 'mag', description = 'Error of the magnification in magnitudes for Lagrange Point 2 parallax') return row_out def main(string = RawCands_table, xx = 0): tt= xxtime() RawCands = Table.read(string, format = 'votable') l = int(len(RawCands)/4) if xx == 0: RawCands = RawCands if xx == 1: RawCands = RawCands[0:l] if xx == 2: RawCands = RawCands[l:2*l] if xx == 3: RawCands = RawCands[2*l:3*l] if xx == 4: RawCands = RawCands[3*l:] table_out = function_to_process(RawCands) table_out.write(Folder+'amlensing' + ['','_a','_b','_c','_d'][xx] +'.vot', format = 'votable',overwrite=True ) print(xxtime()-tt) if xx == 1: subprocess.Popen('python amlensing3.py -t 3', shell=True) if xx == 2: subprocess.Popen('python amlensing3.py -t 4', shell=True) if xx != 0: if os.path.isfile(Folder+'amlensing_a.vot') and os.path.isfile(Folder+'amlensing_b.vot') and os.path.isfile(Folder+'amlensing_c.vot') and os.path.isfile(Folder+'amlensing_d.vot'): processed_sublist = [] for i in range(4): if i +1 == xx: processed_sublist.append(table_out) else: processed_sublist.append(Table.read(Folder+'amlensing'+ ['_a','_b','_c','_d'][i]+'.vot', format = 'votable')) table_out = vstack(processed_sublist) table_out.write(Folder+'amlensing.vot', format = 'votable',overwrite=True) for j in['_a','_b','_c','_d']: os.remove(Folder+'amlensing'+j+'.vot') def check(row): if any(good_hpms['hpms_source_id'] == row['source_id']): if row['ob_ra_error']*row['ob_ra_error'] + row['ob_dec_error']*row['ob_dec_error'] < 100: if (row['pmra']-row['ob_pmra'])*(row['pmra']-row['ob_pmra'])+ (row['pmdec']-row['ob_pmdec'])*(row['pmdec']-row['ob_pmdec']) > 0.7*0.7* (row['pmdec']*row['pmdec']+row['pmra']*row['pmra']): if row['parallax'] > row['ob_parallax']: if row['ob_parallax'] > - 3 * row['ob_parallax_error'] + zeropoint: return True if not row['ob_pmra']: return True return False def function_to_process(sublist): tt=xxtime() table_out = None kk = 0 for row in sublist: kk+=1 if kk%1000 == 0: print(kk,xxtime()-tt) if check(row): row_out = iterrow(row) if row_out is not None: if table_out is None: table_out = row_out else: table_out.add_row(row_out[0]) return table_out if __name__ == '__main__': arg=sys.argv[1:] if len(arg)> 0: if '-t' in ''.join(arg): a = np.where(np.array(arg) == '-t')[0][0]+1 sub = int(arg[a]) good_hpms = Table.read(good_hpms_file, format = 'votable') main(xx = sub ) else: subprocess.Popen('python amlensing3.py -t 1', shell=True) subprocess.Popen('python amlensing3.py -t 2', shell=True)