""" A grammar that turns plc.rawcands to plc.data. This is the "legacy" version used for Svea's 2010 data (in q.rd). Newer data uses makecands.py. """ from __future__ import with_statement import math import os from gavo import api from gavo.base import coords from gavo.grammars import customgrammar def cosdeg(arg): return math.cos(arg/180*math.pi) def sindeg(arg): return math.sin(arg/180*math.pi) def findMinimum(f, left, right, minInterval=3e-8): r = 0.61803399 c = 1-r # set initial parameters x0 = left x3 = right x1 = c*x3 + r*x0 x2 = r*x3 + c*x0 f1 = f(x1) f2 = f(x2) while math.fabs(x3-x0) > minInterval*(math.fabs(x1)+math.fabs(x2)): 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 findMinimum(func, right, left, minInterval=3e-8): # return utils.findMinimum(func, right, left, minInterval) def estimateClosest(lensalpha, lensdelta, obalpha, obdelta, lenspmalpha, lenspmdelta, obpmalpha, obpmdelta, lenspmalphaErr, lenspmdeltaErr): """returns date and distance at the estimated closest approach of lens and object. """ # doing this analytically is a huge PITA AFAICT (minimize the scalar product # of the cartesian motions in the tangential planes -- talk about exploding # terms. Plus, you'll get a grade 4 equation...) dist = lambda t: coords.getGCDist( coords.movePm(lensalpha, lensdelta, lenspmalpha, lenspmdelta, t), coords.movePm(obalpha, obdelta, obpmalpha, obpmdelta, t)) closestTime = findMinimum(dist, -50, 100) closestDist = dist(closestTime) closestDate = 2000+closestTime return closestDate, closestDist POSERR_DEFAULT = 1e-4 PMERR_DEFAULT = 1e-5 def estimateErrors(row, minDate, minDist): DeltaAlpha = (row["alpha"]-row["lensedra"]) DeltaDelta = (row["delta"]-row["lensedde"]) DeltaPmAlpha = row["pma"]/cosdeg(row["delta"] )-row["lensedpmra"]/cosdeg(row["lensedde"]) DeltaPmDelta = row["pmd"]-row["lensedpmde"] cd = cosdeg(row["delta"]) A = DeltaAlpha+DeltaPmAlpha*(minDate-2000) D = DeltaDelta+DeltaPmDelta*(minDate-2000) lensalphaErr = row["alphaErr"] or POSERR_DEFAULT lensdeltaErr = row["deltaErr"] or POSERR_DEFAULT obalphaErr = row["e_lensedra"] obdeltaErr = row["e_lensedde"] lenspmalphaErr = row["pmaErr"] or PMERR_DEFAULT lenspmdeltaErr = row["pmdErr"] or PMERR_DEFAULT obpmalphaErr = row["e_lensedpmra"] obpmdeltaErr = row["e_lensedpmde"] lensdelta = row["delta"] minDateErr = math.sqrt( (pow(cd,4)*pow(DeltaPmAlpha,2)*( pow(lensalphaErr,2)+pow(obalphaErr,2)) +pow(DeltaPmDelta,2)*( pow(lensdeltaErr,2)+pow(obdeltaErr,2))) /pow( pow(cd,2)*pow(DeltaPmAlpha,2)+pow(DeltaPmDelta,2),2) +(pow( pow(cd,2)*DeltaPmAlpha*( DeltaDelta*DeltaPmAlpha -2*DeltaAlpha*DeltaPmDelta) -DeltaDelta*pow(DeltaPmDelta,2),2) *(pow(lenspmdeltaErr,2)+pow(obpmdeltaErr,2)) +pow(cd,4)*pow( DeltaAlpha*( pow(cd,2)*pow(DeltaPmAlpha,2) -pow(DeltaPmDelta,2)) +2*DeltaDelta*DeltaPmAlpha*DeltaPmDelta,2) *(pow(lenspmalphaErr,2) +pow(obpmalphaErr,2))) /pow(pow(cd,2)*pow(DeltaPmAlpha,2)+DeltaPmDelta**2,4)) minDistErr = math.sqrt((pow(A,2)*pow(cd,4)*(pow(lensalphaErr,2)+pow(obalphaErr,2)) +pow(-D-math.pi/90*pow(A,2)*cd*sindeg(lensdelta),2) *pow(lensdeltaErr,2) +pow(D*obdeltaErr,2)+pow((minDate-2000)*A,2)*pow(cd,4)*(pow(lenspmalphaErr,2)+pow( obpmalphaErr,2))+pow((minDate-2000)*D,2)*(pow(lenspmdeltaErr,2)+pow(obpmdeltaErr,2)) +pow(-DeltaPmDelta*D-DeltaPmAlpha*A*pow(cd,2),2)*pow(minDateErr,2)) /(pow(D,2)+pow(A*cd,2))) return minDistErr, minDateErr def makeDataPack(grammar): # the file of confirmed object encounters gcpath = os.path.join(grammar.rd.resdir, "data", "goodCandidates.txt") with open(gcpath) as f: return set( l.strip() for l in f.readlines() if not l.startswith("#")) class RowIterator(customgrammar.CustomRowIterator): def _iterRows(self): with api.getTableConn() as conn: rawcandTD = self.grammar.rd.getById("rawcands") rawcands = api.TableForDef(rawcandTD, connection=conn) for row in rawcands.iterQuery(rawcandTD, None): srcId = row["source"]+" "+row["sourceId"] lensalpha, lensdelta = row["alpha"], row["delta"] lenspmalpha, lenspmdelta = row["pma"], row["pmd"] lensalphaErr, lensdeltaErr = row["alphaErr"], row["deltaErr"] lenspmalphaErr, lenspmdeltaErr = row["pmaErr"], row["pmdErr"] lensmag = row["mag"] lensRmag, lensBmag, lensVmag = row["Rmag"], row["Bmag"], row["Vmag"] lensJmag, lensHmag, lensKmag = row["Jmag"], row["Hmag"], row["Kmag"] obalpha, obdelta = row["lensedra"], row["lensedde"] obpmalpha, obpmdelta = row["lensedpmra"], row["lensedpmde"] obpmalphaErr, obpmdeltaErr = row["e_lensedpmra"], row["e_lensedpmde"] obalphaErr, obdeltaErr = row["e_lensedra"], row["e_lensedde"] #if (lenspmalpha+obpmalpha)>1e-10 and (lenspmdelta+obpmdelta)>1e-10: #if abs(lenspmalpha-obpmalpha)/abs(lenspmalpha+obpmalpha)<0.05 and abs(lenspmdelta-obpmdelta)/abs(lenspmdelta+obpmdelta)<0.05: #continue minDate, minDist = estimateClosest(lensalpha, lensdelta, obalpha, obdelta, lenspmalpha, lenspmdelta, obpmalpha, obpmdelta, lenspmalphaErr, lenspmdeltaErr) minDistErr, minDateErr = estimateErrors(row, minDate, minDist) obJ = row["lensedjmag"] obH = row["lensedhmag"] obK = row["lensedkmag"] obB = row["lensedBmag"] obR = row["lensedRmag"] obI = row["lensedImag"] confirmed = srcId in self.grammar.dataPack source = row["source"] yield locals() # vi:et:sta:sw=2:ts=2:ai