import numpy as np import struct import pylab import math import os import time import random from scipy import stats from matplotlib.backends.backend_pdf import FigureCanvasPdf as FigureCanvas from matplotlib.figure import Figure def anomcalc(meananom,ecc): """calculates eccentric anomaly from the mean anomaly""" #based on Murray and Dermot meananom=np.asarray(meananom) ecc=np.asarray(ecc) k=0.85 sign=np.sin(meananom)/np.fabs(np.sin(meananom)) E0=meananom+sign*k*ecc E=E0 sel=np.arange(0,len(E0)) fracdiff=sel*0.+1. sel2=np.where(ecc >= 1) fracdiff[sel2]=0. while np.amax(fracdiff) > 0.001: f0=E[sel]-ecc[sel]*np.sin(E[sel])-meananom[sel] f1=1.-ecc[sel]*np.cos(E[sel]) f2=ecc[sel]*np.sin(E[sel]) f3=ecc[sel]*np.cos(E[sel]) del1=-f0/f1 del2=-f0/(f1+.5*del1*f2) del3=-f0/(f1+.5*del2*f2+1./6.*del2**2*f3) nextone=E[sel]+del3 diff=E[sel]-nextone fracdiff[sel]=np.fabs(diff/E[sel]) E[sel]=nextone sel=np.where((fracdiff > 0.001)) sel2=np.where(ecc >= 1) fracdiff[sel2]=0. sel2=np.where(ecc >= 1) E[sel2]=2.*math.pi*1.0001 return(E) def angcalc(E,ecc): """calculates the true anomaly from the eccentric anomaly""" E=np.asarray(E) ecc=np.asarray(ecc) cosang=(np.cos(E)-ecc)/(1-ecc*np.cos(E)) angpos=np.arccos(cosang) #sometimes E is just under 0...this is wrong E=np.fabs(E) #output of arccos goes from 0 to pi, but real angles go up to 2*pi #we need to correct for this acoscorrec=np.floor(E/math.pi)*2*math.pi #however sometime E is just over 2*pi. these don't need corrected Enoisecorrec=np.floor(E/(2*math.pi))*4.*math.pi acoscorrec=acoscorrec-Enoisecorrec #correcting angpos angpos=acoscorrec-angpos angpos=(angpos*angpos)**.5 return(angpos) def rcalc(E,a,ecc): """calculates distance from Sun for a given orbit and ecc. anomaly""" r=a*(1.-ecc*np.cos(E)) return(r) def orb2xv(a,e,inc,capom,omega,capm,mcent): """converts from orbital elements to cartesian""" #transform mean anomalies into true anomaly and distance eccanom=anomcalc(capm,e) angpos=angcalc(eccanom,e) rpos=rcalc(eccanom,a,e) #transform into actual positions and velocities #build transformation matrix transformx = (np.cos(capom) * np.cos(omega + angpos) - np.sin(capom) * np.sin(omega + angpos) * np.cos(inc)) transformy = (np.sin(capom) * np.cos(omega + angpos) + np.cos(capom) * np.sin(omega + angpos) * np.cos(inc)) transformz = np.sin(omega + angpos) * np.sin(inc) #transforming coords xg = rpos * transformx yg = rpos * transformy zg = rpos * transformz transformx, transformy, transformz = None, None, None #calculating velocities in specific orbit coordinate period = np.sqrt(4 * math.pi**2 * a**3 / mcent) angvel = 2 * math.pi / period #doing velocity transformations...see p. 31 and 51 of Murray and Dermott xdot = -angvel * a * np.sin(angpos) / np.sqrt(1 - e**2) ydot = angvel * a * (e + np.cos(angpos)) / np.sqrt(1 - e**2) zdot = 0 * angpos vxg = ((np.cos(omega) * np.cos(capom) - np.sin(omega) * np.sin(capom) * np.cos(inc)) * xdot - (np.sin(omega) * np.cos(capom) + np.cos(omega) * np.sin(capom) * np.cos(inc)) * ydot) vyg = ((np.cos(omega) * np.sin(capom) + np.sin(omega) * np.cos(capom) * np.cos(inc)) * xdot + (np.cos(omega) * np.cos(capom) * np.cos(inc) - np.sin(omega) * np.sin(capom)) * ydot) vzg = ((np.sin(omega) * np.sin(inc)) * xdot + (np.cos(omega) * np.sin(inc)) * ydot) xdot, ydot, zdot = None, None, None rpos, angpos, angvel, period, eccanom = None, None, None, None, None return(xg,yg,zg,vxg,vyg,vzg) def xv2orb(x,y,z,vx,vy,vz,mcent): """this program converts xv coordinates to orbital elements""" tiny=4e-15 dr = 180.0/math.pi #calculate angular momenta hx = y * vz - z * vy hy = z * vx - x * vz hz = x * vy - y * vx h2 = hx * hx + hy * hy + hz * hz h = np.sqrt(h2) inc = np.arccos(hz / h) u = np.empty(len(x)) capom = np.empty(len(x)) omega = np.empty(len(x)) capm = np.empty(len(x)) e = np.empty(len(x)) a = np.empty(len(x)) #calculating long. of ascending node fac = np.sqrt(hx**2 + hy**2) / h sel = np.unique(np.asarray(np.where(fac < tiny))) capom[sel] = 0.0 u[sel] = np.arctan2(y[sel], x[sel]) sel2 = np.unique(np.asarray(np.where(np.fabs(inc[sel] - math.pi) < 10.0 * tiny))) u[sel[sel2]] = -u[sel[sel2]] sel = np.unique(np.asarray(np.where(fac >= tiny))) capom[sel] = np.arctan2(hx[sel], -hy[sel]) u[sel] = (np.arctan2((z[sel] / np.sin(inc[sel])), (x[sel] * np.cos(capom[sel]) + y[sel] * np.sin(capom[sel])))) sel = np.unique(np.asarray(np.where(capom < 0))) capom[sel] = capom[sel] + 2 * math.pi sel = np.unique(np.asarray(np.where(u < 0))) u[sel] = u[sel] + 2 * math.pi #calculating energy, etc. r = np.sqrt(x * x + y * y + z * z) v2 = vx * vx + vy * vy + vz * vz v = np.sqrt(v2) vdotr = x * vx + y * vy + z * vz energy = v2 / 2 - mcent / r #determing conic section and label it with ialpha ialpha = np.empty(len(x), dtype = int) sel = np.unique(np.asarray(np.where(energy < 0))) ialpha[sel] = -1 sel = np.unique(np.asarray(np.where(energy > 0))) ialpha[sel] = 1 sel = np.unique(np.asarray(np.where(np.fabs(energy * r / mcent) < np.sqrt(tiny)))) ialpha[sel] = 0 #ellipse ellipse = np.unique(np.asarray(np.where(ialpha == -1))) a[ellipse] = -mcent / energy[ellipse] / 2 fac = 1 - h2[ellipse] / mcent / a[ellipse] w = np.empty(len(fac)) cape = np.empty(len(fac)) sel = np.unique(np.asarray(np.where(fac > tiny))) e[ellipse[sel]] = np.sqrt(fac[sel]) face = ((a[ellipse[sel]] - r[ellipse[sel]]) / (a[ellipse[sel]] * e[ellipse[sel]])) sel2 = np.unique(np.asarray(np.where(face > 1))) cape[sel[sel2]] = 0.0 sel2 = np.unique(np.asarray(np.where((face > -1) & (face <= 1)))) cape[sel[sel2]] = np.arccos(face[sel2]) sel2 = np.unique(np.asarray(np.where(face <= -1))) cape[sel[sel2]] = math.pi sel2 = np.unique(np.asarray(np.where(vdotr[ellipse[sel]] < 0))) cape[sel[sel2]] = 2 * math.pi - cape[sel[sel2]] cw = np.empty(len(sel)) sw = np.empty(len(sel)) cw = ((np.cos(cape[sel]) - e[ellipse[sel]])/ (1 - e[ellipse[sel]] * np.cos(cape[sel]))) sw = (np.sqrt(1 - e[ellipse[sel]]**2) * np.sin(cape[sel]) / (1 - e[ellipse[sel]] * np.cos(cape[sel]))) w[sel] = np.arctan2(sw, cw) sel2 = np.unique(np.asarray(np.where(w[sel] < 0))) w[sel[sel2]] = w[sel[sel2]] + 2 * math.pi #taking care of cases of almost perfectly circular orbits sel = np.unique(np.asarray(np.where(fac <= tiny))) e[ellipse[sel]] = 0.0 w[sel] = u[ellipse[sel]] cape[sel] = u[ellipse[sel]] capm[ellipse] = cape - e[ellipse] * np.sin(cape) omega[ellipse] = u[ellipse] - w sel = np.unique(np.asarray(np.where(omega[ellipse] < 0))) omega[ellipse[sel]] = omega[ellipse[sel]] + 2 * math.pi omega[ellipse] = (omega[ellipse] - np.floor(omega[ellipse] / 2 / math.pi) * 2 * math.pi) #hyperbola hyper = np.unique(np.asarray(np.where(ialpha == 1))) a[hyper] = mcent / energy[hyper] / 2 fac = h2[hyper] / (mcent * a[hyper]) w = np.empty(len(fac)) tmpf = np.empty(len(fac)) capf = np.empty(len(fac)) sel = np.unique(np.asarray(np.where(fac > tiny))) e[hyper[sel]] = np.sqrt(1 + fac[sel]) tmpf[sel] = ((a[hyper[sel]] + r[hyper[sel]]) / (a[hyper[sel]] * e[hyper[sel]])) sel2 = np.unique(np.asarray(np.where(tmpf[sel] < 1))) tmpf[sel[sel2]] = 1 capf[sel] = np.log10(tmpf[sel] + np.sqrt(tmpf[sel]**2 - 1)) sel2 = np.unique(np.asarray(np.where(vdotr[hyper[sel]] < 0))) capf[sel[sel2]] = -capf[sel[sel2]] cw = ((e[hyper[sel]] - np.cosh(capf[sel])) / (e[hyper[sel]] * np.cosh(capf[sel]))) sw = ((e[hyper[sel]]**2 - 1) * np.sinh(capf[sel]) / (e[hyper[sel]] * np.cosh(capf[sel]))) w[sel] = np.arctan2(sw, cw) sel2 = np.unique(np.asarray(np.where(w[sel] < 0))) w[sel[sel2]] = w[sel[sel2]] + 2 * math.pi #taking care of near-parabola cases sel = np.unique(np.asarray(np.where(fac <= tiny))) e[hyper[sel]] = 1 tmpf[sel] = h2[hyper[sel]] / mcent / 2 w[sel] = np.arccos(2 * tmpf[sel] / r[hyper[sel]] -1) sel2 = np.unique(np.asarray(np.where(vdotr[hyper[sel]] < 0))) w[sel[sel2]] = 2 * math.pi - w[sel[sel2]] tmpf[sel] = ((a[hyper[sel]] + r[hyper[sel]]) / (a[hyper[sel]] * e[hyper[sel]])) capf[sel] = np.log10(tmpf[sel] + np.sqrt(tmpf[sel]**2 - 1)) capm[hyper] = e[hyper] * np.sinh(capf) - capf omega[hyper] = u[hyper] - w sel = np.unique(np.asarray(np.where(omega[hyper] < 0))) omega[hyper[sel]] = omega[hyper[sel]] + 2 * math.pi omega[hyper] = (omega[hyper] - np.floor(omega[hyper] / 2 / math.pi) * 2 * math.pi) #parabola para = np.unique(np.asarray(np.where(ialpha == 0))) a[para] = h2[para] / mcent / 2 e[para] = 1 w = np.arccos(2 * a[para] / r[para] - 1) sel = np.unique(np.asarray(np.where(vdotr[para] < 0))) w[sel] = 2 * math.pi - w[sel] tmpf = np.tan(w / 2) capm[para] = tmpf * (1 + tmpf**2 / 3) omega[para] = u[para] - w sel = np.unique(np.asarray(np.where(omega[para] < 0))) omega[para[sel]] = omega[para[sel]] + 2 * math.pi omega[para] = (omega[para] - np.floor(omega[para] / 2 / math.pi) * 2 * math.pi) return(a,e,inc,capom,omega,capm) def anomcalc(meananom,ecc): """calculates eccentric anomaly from the mean anomaly""" #based on Murray and Dermot meananom=np.asarray(meananom) ecc=np.asarray(ecc) k=0.85 sign=np.sin(meananom)/np.fabs(np.sin(meananom)) E0=meananom+sign*k*ecc E=E0 sel=np.arange(0,len(E0)) fracdiff=sel*0.+1. sel2=np.where(ecc >= 1) fracdiff[sel2]=0. while np.amax(fracdiff) > 0.001: f0=E[sel]-ecc[sel]*np.sin(E[sel])-meananom[sel] f1=1.-ecc[sel]*np.cos(E[sel]) f2=ecc[sel]*np.sin(E[sel]) f3=ecc[sel]*np.cos(E[sel]) del1=-f0/f1 del2=-f0/(f1+.5*del1*f2) del3=-f0/(f1+.5*del2*f2+1./6.*del2**2*f3) nextone=E[sel]+del3 diff=E[sel]-nextone fracdiff[sel]=np.fabs(diff/E[sel]) E[sel]=nextone sel=np.where((fracdiff > 0.001)) sel2=np.where(ecc >= 1) fracdiff[sel2]=0. sel2=np.where(ecc >= 1) E[sel2]=2.*math.pi*1.0001 return(E) def angcalc(E,ecc): """calculates the true anomaly from the eccentric anomaly""" E=np.asarray(E) ecc=np.asarray(ecc) cosang=(np.cos(E)-ecc)/(1-ecc*np.cos(E)) angpos=np.arccos(cosang) #sometimes E is just under 0...this is wrong E=np.fabs(E) #output of arccos goes from 0 to pi, but real angles go up to 2*pi #we need to correct for this acoscorrec=np.floor(E/math.pi)*2*math.pi #however sometime E is just over 2*pi. these don't need corrected Enoisecorrec=np.floor(E/(2*math.pi))*4.*math.pi acoscorrec=acoscorrec-Enoisecorrec #correcting angpos angpos=acoscorrec-angpos angpos=(angpos*angpos)**.5 return(angpos) def rcalc(E,a,ecc): """calculates distance from Sun for a given orbit and ecc. anomaly""" r=a*(1.-ecc*np.cos(E)) return(r) def chunks1(string,n,offset,width): return [string[i+offset:i+offset+width] for i in range(0,len(string),n)] def chunks2(string,n,linesize,offset,width): return [string[i+offset+j*linesize:i+offset+width+j*linesize] for i in range(0,len(string),n) for j in range(4)] idall = np.empty(0, dtype=int) tall = np.empty(0) aall = np.empty(0) incall = np.empty(0) omegaall = np.empty(0) rdotvall = np.empty(0) incecall = np.empty(0) omegaecall = np.empty(0) capomecall = np.empty(0) capomall = np.empty(0) qall = np.empty(0) oldk=-1 fpos = 0 k = 0 ndir = 1000 print('reading flux files') while (k < ndir): filename = str(k+1)+'/flux.dat' while ((os.path.exists(filename)==False)&(koldk): f = open(filename, 'rb') linesize = 38 entrysize = 16+linesize chunk = entrysize*1000000 bytesread = min(fsize-fpos,chunk) a = f.read(bytesread) #compile times headertimest = chunks1(a,entrysize,4,4) time = np.asarray([struct.unpack('f',timestampst)[0] for timestampst in headertimest]) #compile tp data idst = chunks1(a,entrysize,20,2) mst = chunks1(a,entrysize,22,4) xst = chunks1(a,entrysize,26,4) yst = chunks1(a,entrysize,30,4) zst = chunks1(a,entrysize,34,4) vxst = chunks1(a,entrysize,38,4) vyst = chunks1(a,entrysize,42,4) vzst = chunks1(a,entrysize,46,4) idtp = np.asarray([struct.unpack('h',idst1)[0] for idst1 in idst]) xtp = np.asarray([struct.unpack('f',xst1)[0] for xst1 in xst]) ytp = np.asarray([struct.unpack('f',yst1)[0] for yst1 in yst]) ztp = np.asarray([struct.unpack('f',zst1)[0] for zst1 in zst]) vxtp = np.asarray([struct.unpack('f',vxst1)[0] for vxst1 in vxst]) vytp = np.asarray([struct.unpack('f',vyst1)[0] for vyst1 in vyst]) vztp = np.asarray([struct.unpack('f',vzst1)[0] for vzst1 in vzst]) #calculate barycenter for each passage nentries = len(idtp) xbar = np.zeros(nentries) ybar = np.zeros(nentries) zbar = np.zeros(nentries) vxbar = np.zeros(nentries) vybar = np.zeros(nentries) vzbar = np.zeros(nentries) mcent = 2.9630927804939294e-4 xb = xtp - xbar yb = ytp - ybar zb = ztp - zbar vxb = vxtp - vxbar vyb = vytp - vybar vzb = vztp - vzbar a, e, inc, capom, omega, capm = xv2orb(xb, yb, zb, vxb, vyb, vzb,mcent) q = a*(1.-e) galang = 60.2 * math.pi/180. newx = xb*np.cos(galang)-zb*np.sin(galang) newy = yb newz = xb*np.sin(galang)+zb*np.cos(galang) newvx = vxb*np.cos(galang)-vzb*np.sin(galang) newvy = vyb newvz = vxb*np.sin(galang)+vzb*np.cos(galang) xbar = None ybar = None zbar = None vxbar = None vybar = None vzbar = None #Rotate to the ecliptic #see if it's comet in or going out rdotv = np.sign(newx*newvx + newy*newvy + newz*newvz) a, e, incec, capomec, omegaec, capm = xv2orb(newx, newy, newz, newvx, newvy, newvz,mcent) q = a*(1.-e) sel = np.where(a>0)[0] idall = np.concatenate((idall,idtp[sel]+k*100)) qall = np.concatenate((qall,q[sel])) aall = np.concatenate((aall,a[sel])) incall = np.concatenate((incall,inc[sel])) omegaall = np.concatenate((omegaall,omega[sel])) capomall = np.concatenate((capomall,capom[sel])) tall = np.concatenate((tall,time[sel])) incecall = np.concatenate((incecall,incec[sel])) omegaecall = np.concatenate((omegaecall,omegaec[sel])) capomecall = np.concatenate((capomecall,capomec[sel])) rdotvall = np.concatenate((rdotvall,rdotv[sel])) fpos = f.tell() oldk = k if (fpos==fsize): f.close() k = k+1 fpos = 0 print('done reading flux files') #only keep LPCs on incoming trajectories incoming = np.where(rdotvall<0)[0] idall = idall[incoming] tall = tall[incoming] aall = aall[incoming] incall = incall[incoming] omegaall = omegaall[incoming] rdotvall = rdotvall[incoming] incecall = incecall[incoming] omegaecall = omegaecall[incoming] capomecall = capomecall[incoming] capomall = capomall[incoming] qall = qall[incoming] #only keep LPCs that enter inner 20 AU after t=500 Myrs #and that make at least one pericenter passage between 8 and 12 AU print('assembling bad id list') iduniq = np.unique(idall) badids = [] nonnearids = [] for idval in iduniq: sel = np.where(idall==idval)[0] #test entry time tvals = np.sort(tall[sel]) if tvals[0]/365e6 < 500: badids.append(idval) nonnearids.append(idval) #search for pericenter passage obsanalogs = np.where((qall[sel]>8)&(qall[sel]<12))[0] if len(obsanalogs)==0: badids.append(idval) nearanalogs = np.where((qall[sel]>0)&(qall[sel]<4))[0] if len(nearanalogs)==0: nonnearids.append(idval) print('done assembling bad id list') #strip out bad id data print('stripping out bad ids') selbadall = np.empty(0,dtype='int').flatten() for badid in badids: selbad = np.where(idall==badid)[0] selbadall = np.concatenate((selbadall,selbad.flatten())) selfarall = np.empty(0,dtype='int').flatten() for badid in nonnearids: selfar = np.where(idall==badid)[0] selfarall = np.concatenate((selfarall,selfar.flatten())) idnear = np.delete(idall,selfarall) qnear = np.delete(qall,selfarall) anear = np.delete(aall,selfarall) incnear = np.delete(incall,selfarall) omeganear = np.delete(omegaall,selfarall) capomnear = np.delete(capomall,selfarall) tnear = np.delete(tall,selfarall) incecnear = np.delete(incecall,selfarall) omegaecnear = np.delete(omegaecall,selfarall) capomecnear = np.delete(capomecall,selfarall) rdotvnear = np.delete(rdotvall,selfarall) idall = np.delete(idall,selbadall) qall = np.delete(qall,selbadall) aall = np.delete(aall,selbadall) incall = np.delete(incall,selbadall) omegaall = np.delete(omegaall,selbadall) capomall = np.delete(capomall,selbadall) tall = np.delete(tall,selbadall) incecall = np.delete(incecall,selbadall) omegaecall = np.delete(omegaecall,selbadall) capomecall = np.delete(capomecall,selbadall) rdotvall = np.delete(rdotvall,selbadall) print('done stripping out bad ids') #make sure everything is ordered properly with time print('sorting') order = np.argsort(tall) idall = idall[order] qall = qall[order] aall = aall[order] incall = incall[order] omegaall = omegaall[order] capomall = capomall[order] tall = tall[order] incecall = incecall[order] omegaecall = omegaecall[order] capomecall = capomecall[order] rdotvall = rdotvall[order] order = np.argsort(tnear) idnear = idnear[order] qnear = qnear[order] anear = anear[order] incnear = incnear[order] omeganear = omeganear[order] capomnear = capomnear[order] tnear = tnear[order] incecnear = incecnear[order] omegaecnear = omegaecnear[order] capomecnear = capomecnear[order] rdotvnear = rdotvnear[order] print('done sorting') #build longitudes of perihelion pomegaecall = np.mod(omegaecall+capomecall,2.*math.pi) pomegaall = np.mod(omegaall+capomall,2.*math.pi) pomegaecnear = np.mod(omegaecnear+capomecnear,2.*math.pi) pomeganear = np.mod(omeganear+capomnear,2.*math.pi) id0,frac,a0,e0,inc0,capom0,omega0,capm0 = np.genfromtxt('populations/parts0.dat',unpack=True,skip_header=1) a0all = np.zeros(len(idall))-1. a0near = np.zeros(len(idnear))-1. #build npass array #first column is passages inside 8.5 #last column is passages inside 20 npassall = np.zeros((24,len(idall))) npassnear = np.zeros((24,len(idnear))) qcritvals = np.arange(8.5,20.5,0.5) #go through each comet and add up passages within a certain pericenter value print('assigning passages') iduniq = np.unique(idall) tlast = np.zeros(len(idall)) for k,idval in enumerate(iduniq): #select all peri passages for a single comet comsel = np.sort(np.where(idall==idval)[0]) #cycle through the list of q vals for i,qval in enumerate(qcritvals): sel = np.sort(np.where(qall[comsel] 0: sel = np.sort(sel.flatten()) for j,indice in enumerate(sel): #increase passages by 1 for current index and all following ones npassall[i,comsel[indice:]] = npassall[i,comsel[indice:]] + 1 #also, assign initial semimajor axes ind = np.where(id0==idval)[0] a0all[comsel] = a0[ind] iduniq = np.unique(idnear) for idval in iduniq: #select all peri passages for a single comet comsel = np.sort(np.where(idnear==idval)[0]) #cycle through the list of q vals for i,qval in enumerate(qcritvals): sel = np.sort(np.where(qnear[comsel] 0: sel = np.sort(sel.flatten()) for j,indice in enumerate(sel): #increase passages by 1 for current index and all following ones npassnear[i,comsel[indice:]] = npassnear[i,comsel[indice:]] + 1 #also, assign initial semimajor axes ind = np.where(id0==idval)[0] a0near[comsel] = a0[ind] print('done assigning passages') print('assigning previous perihelion passages') qprevall = np.zeros(len(idall)) wrapperlen = 4 linelen = 34 + wrapperlen * 2 headerlen = 8 + wrapperlen * 2 for i,idval in enumerate(idall): #build path name and open binary file path = 'backintegrationsallapps/' + str(idall[i]) + '/' + str(int(tall[i])) filename = path+'/bin.dat' f = open(filename, 'rb') f.seek(linelen + headerlen*2,0) a = f.read(linelen) #output binary info ast = a[wrapperlen + 10:wrapperlen + 14] est = a[wrapperlen + 14:wrapperlen + 18] #unpack binary data aprev = np.asarray(struct.unpack('f',ast)) eprev = np.asarray(struct.unpack('f',est)) #calculate previous perihelion qprevall[i] = aprev * (1. - eprev) print('done assigning previous perihelion passages') print('reading observed LPC data') exec(open('/Users/nathankaib/Work/RealLPCs/converttogalcoords.py').read()) print('done reading observed LPC data') print('reading observed LPC back integrations') #assemble list of LPC names comname = np.genfromtxt('/Users/nathankaib/Work/RealLPCs/distantLPCs_JPL.txt',unpack=True,usecols=(0),skip_header=1,delimiter=',',dtype='str') for i,name in enumerate(comname): ind1 = name.find('C') ind2 = name.find('(') name = name[ind1:ind2-1] name = name.replace('/','') name = name.replace(' ','_') comname[i] = name #read in the binary data for each comet back integration qprevreal = np.zeros(len(comname)) wrapperlen = 4 linelen = 34 + wrapperlen * 2 headerlen = 8 + wrapperlen * 2 for i,comnameval in enumerate(comname): #build path name and open binary file path = '/Users/nathankaib/Work/RealLPCs/backintegrations/'+comname[i] filename = path+'/bin.dat' f = open(filename, 'rb') f.seek(linelen + headerlen*2,0) a = f.read(linelen) #output binary info ast = a[wrapperlen + 10:wrapperlen + 14] est = a[wrapperlen + 14:wrapperlen + 18] #unpack binary data aprev = np.asarray(struct.unpack('f',ast)) eprev = np.asarray(struct.unpack('f',est)) #calculate previous perihelion qprevreal[i] = aprev * (1. - eprev) print('done reading observed LPC back integrations') #now go through and find the min and max 2-sigma Npass value at each qcrit value npassvals = np.arange(1,20.5,0.5) apvals = np.zeros((len(npassvals),len(qcritvals))) omegapvals = np.zeros((len(npassvals),len(qcritvals))) qprevpvals = np.zeros((len(npassvals),len(qcritvals))) apvals_acap = np.zeros((len(npassvals),len(qcritvals))) omegapvals_acap = np.zeros((len(npassvals),len(qcritvals))) qprevpvals_acap = np.zeros((len(npassvals),len(qcritvals))) for i,qval in enumerate(qcritvals): print(i) for j,npassval in enumerate(npassvals): sel = np.where((qall>8)&(qall<12)&(npassall[i,:]<=npassval))[0] #account for cases where npassval is not a whole number selbubble = np.where((qall>8)&(qall<12)&(npassall[i,:]npassval))[0] if len(selbubble)>0: deciders = np.random.random_sample(len(selbubble)) probcut = npassval - int(npassval) keep = np.where(deciders<=probcut)[0] if len(keep)>0: sel = np.concatenate((sel.flatten(),selbubble[keep].flatten())) #use K-S test to assess fitness if len(sel)>0: d,apvals[j,i] = stats.ks_2samp(ag,aall[sel]) d,omegapvals[j,i] = stats.ks_2samp(np.sin(2.*omegag*math.pi/180.),np.sin(omegaall[sel]*2.)) d,qprevpvals[j,i] = stats.ks_2samp(qprevreal,qprevall[sel]) sel = np.where((qall>8)&(qall<12)&(npassall[i,:]<=npassval)&(a0all<2.5e4))[0] #account for cases where npassval is not a whole number selbubble = np.where((qall>8)&(qall<12)&(npassall[i,:]npassval)&(a0all<2.5e4))[0] if len(selbubble)>0: deciders = np.random.random_sample(len(selbubble)) probcut = npassval - int(npassval) keep = np.where(deciders<=probcut)[0] if len(keep)>0: sel = np.concatenate((sel.flatten(),selbubble[keep].flatten())) #use K-S test to assess fitness if len(sel)>0: d,apvals_acap[j,i] = stats.ks_2samp(ag,aall[sel]) d,omegapvals_acap[j,i] = stats.ks_2samp(np.sin(2.*omegag*math.pi/180.),np.sin(omegaall[sel]*2.)) d,qprevpvals_acap[j,i] = stats.ks_2samp(qprevreal,qprevall[sel]) figpath = 'SMAMatch.pdf' fig = Figure(figsize=(8,6)) canvas = FigureCanvas(fig) logbins = np.logspace(2.5,5.01,10000) ax = fig.add_subplot(1,1,1) ax.hist(ag[0:],bins=logbins,histtype='step',linewidth=5,alpha=0.5,color='k',range=[1,1e6],density=True,cumulative=True) line1, = ax.plot([0,1],[2,2],linewidth=5,alpha=0.5,color='k') sel = np.where((qall>8)&(qall<12)&(npassall[23,:]<2e5))[0] ax.hist(aall[sel],bins=logbins,histtype='step',linewidth=1,linestyle='dotted',color='k',range=[1,1e6],density=True,cumulative=True) line2, = ax.plot([0,1],[2,2],linewidth=1,linestyle='dotted',color='k') sel = np.where((qall>8)&(qall<12)&(npassall[23,:]<2))[0] ax.hist(aall[sel],bins=logbins,histtype='step',linewidth=1,linestyle='dashed',color='k',range=[1,1e6],density=True,cumulative=True) line3, = ax.plot([0,1],[2,2],linewidth=1,linestyle='dashed',color='k') sel = np.where((qall>8)&(qall<12)&(npassall[23,:]<6))[0] ax.hist(aall[sel],bins=logbins,histtype='step',linewidth=1,linestyle='dashdot',color='k',range=[1,1e6],density=True,cumulative=True) line4, = ax.plot([0,1],[2,2],linewidth=1,linestyle='dashdot',color='k') #sel = np.where((qall>8)&(qall<12)&(npassall>0)&(npassall<6)&(a0all<3e4))[0] #ax.hist(aall[sel],bins=logbins,histtype='step',linewidth=1,linestyle='dashdot',color='k',range=[1,1e6],density=True,cumulative=True) #line5, = ax.plot([0,1],[2,2],linewidth=1,linestyle='dashdot',color='k') #lineall = [line1, line2, line3, line4,line5] #laball = ['Observed',r'N$_{\rm pass}$ = $\infty$',r'N$_{\rm pass}$ = 1',r'N$_{\rm pass}$ = 5',r'N$_{\rm pass}$ = 5, a$_0$ < 30,000 AU'] lineall = [line1, line2, line3, line4] laball = ['Observed',r'N$_{\rm pass}$ = $\infty$',r'N$_{\rm pass}$ = 1',r'N$_{\rm pass}$ = 5'] ax.set_xlabel('a (au)') ax.set_xscale('log') ax.set_ylabel('Cumulative Fraction') ax.set_xlim(4e2,1e5) ax.set_ylim(0,1) ax.legend(lineall,laball,loc='upper left') fig.tight_layout() canvas.print_figure(figpath) figpath = 'SMAMatch_talk1.pdf' fig = Figure(figsize=(8,6)) canvas = FigureCanvas(fig) logbins = np.logspace(2.5,5.01,10000) ax = fig.add_subplot(1,1,1) ax.hist(ag[0:],bins=logbins,histtype='step',linewidth=5,alpha=0.5,color='k',range=[1,1e6],density=True,cumulative=True) line1, = ax.plot([0,1],[2,2],linewidth=5,alpha=0.5,color='k') sel = np.where((qall>8)&(qall<12)&(npassall[23,:]<2e5))[0] ax.hist(aall[sel],bins=logbins,histtype='step',linewidth=1,linestyle='dotted',color='r',range=[1,1e6],density=True,cumulative=True) line2, = ax.plot([0,1],[2,2],linewidth=1,linestyle='dotted',color='r') #sel = np.where((qall>8)&(qall<12)&(npassall[23,:]<2))[0] #ax.hist(aall[sel],bins=logbins,histtype='step',linewidth=1,linestyle='dashed',color='k',range=[1,1e6],density=True,cumulative=True) #line3, = ax.plot([0,1],[2,2],linewidth=1,linestyle='dashed',color='k') #sel = np.where((qall>8)&(qall<12)&(npassall[23,:]<6))[0] #ax.hist(aall[sel],bins=logbins,histtype='step',linewidth=1,linestyle='dashdot',color='k',range=[1,1e6],density=True,cumulative=True) #line4, = ax.plot([0,1],[2,2],linewidth=1,linestyle='dashdot',color='k') #sel = np.where((qall>8)&(qall<12)&(npassall>0)&(npassall<6)&(a0all<3e4))[0] #ax.hist(aall[sel],bins=logbins,histtype='step',linewidth=1,linestyle='dashdot',color='k',range=[1,1e6],density=True,cumulative=True) #line5, = ax.plot([0,1],[2,2],linewidth=1,linestyle='dashdot',color='k') #lineall = [line1, line2, line3, line4,line5] #laball = ['Observed',r'N$_{\rm pass}$ = $\infty$',r'N$_{\rm pass}$ = 1',r'N$_{\rm pass}$ = 5',r'N$_{\rm pass}$ = 5, a$_0$ < 30,000 AU'] lineall = [line1, line2] laball = ['Observed',r'N$_{\rm pass}$ = $\infty$'] ax.set_xlabel('a (au)') ax.set_xscale('log') ax.set_ylabel('Cumulative Fraction') ax.set_xlim(4e2,1e5) ax.set_ylim(0,1) ax.legend(lineall,laball,loc='upper left') figpath = 'SMAMatch_talk2.pdf' fig = Figure(figsize=(8,6)) canvas = FigureCanvas(fig) logbins = np.logspace(2.5,5.01,10000) ax = fig.add_subplot(1,1,1) ax.hist(ag[0:],bins=logbins,histtype='step',linewidth=5,alpha=0.5,color='k',range=[1,1e6],density=True,cumulative=True) line1, = ax.plot([0,1],[2,2],linewidth=5,alpha=0.5,color='k') sel = np.where((qall>8)&(qall<12)&(npassall[23,:]<2e5))[0] ax.hist(aall[sel],bins=logbins,histtype='step',linewidth=1,linestyle='dotted',color='r',range=[1,1e6],density=True,cumulative=True) line2, = ax.plot([0,1],[2,2],linewidth=1,linestyle='dotted',color='r') #sel = np.where((qall>8)&(qall<12)&(npassall[23,:]<2))[0] #ax.hist(aall[sel],bins=logbins,histtype='step',linewidth=1,linestyle='dashed',color='k',range=[1,1e6],density=True,cumulative=True) #line3, = ax.plot([0,1],[2,2],linewidth=1,linestyle='dashed',color='k') sel = np.where((qall>8)&(qall<12)&(npassall[23,:]<6))[0] ax.hist(aall[sel],bins=logbins,histtype='step',linewidth=1,linestyle='dashdot',color='b',range=[1,1e6],density=True,cumulative=True) line4, = ax.plot([0,1],[2,2],linewidth=1,linestyle='dashdot',color='b') #sel = np.where((qall>8)&(qall<12)&(npassall>0)&(npassall<6)&(a0all<3e4))[0] #ax.hist(aall[sel],bins=logbins,histtype='step',linewidth=1,linestyle='dashdot',color='k',range=[1,1e6],density=True,cumulative=True) #line5, = ax.plot([0,1],[2,2],linewidth=1,linestyle='dashdot',color='k') lineall = [line1, line2, line4] laball = ['Observed',r'N$_{\rm pass}$ = $\infty$',r'N$_{\rm pass}$ = 5'] ax.set_xlabel('a (au)') ax.set_xscale('log') ax.set_ylabel('Cumulative Fraction') ax.set_xlim(4e2,1e5) ax.set_ylim(0,1) ax.legend(lineall,laball,loc='upper left') fig.tight_layout() canvas.print_figure(figpath) fig.tight_layout() canvas.print_figure(figpath) figpath = 'OmegaMatch.pdf' fig = Figure(figsize=(8,12)) canvas = FigureCanvas(fig) ax = fig.add_subplot(2,1,1) ax.hist(np.sin(omegag[0:]*2.*math.pi/180.),bins=10000,histtype='step',linewidth=5,color='k',range=[-1,1.01],density=True,cumulative=True,alpha=0.5) line1, = ax.plot([0,1],[2,2],linewidth=5,color='k',alpha=0.5) sel = np.where((qall>8)&(qall<12)&(npassall[23,:]<2e5))[0] ax.hist(np.sin(2.*omegaall[sel]),bins=10000,histtype='step',linewidth=1,linestyle='dotted',color='k',range=[-1,1.01],density=True,cumulative=True) line2, = ax.plot([0,1],[2,2],linewidth=1,linestyle='dotted',color='k') sel = np.where((qall>8)&(qall<12)&(npassall[23,:]<6))[0] ax.hist(np.sin(2.*omegaall[sel]),bins=10000,histtype='step',linewidth=1,linestyle='dashed',color='k',range=[-1,1.01],density=True,cumulative=True) line3, = ax.plot([0,1],[2,2],linewidth=1,linestyle='dashed',color='k') sel = np.where((qall>8)&(qall<12)&(npassall[23,:]<6)&(a0all<4e4))[0] ax.hist(np.sin(2.*omegaall[sel]),bins=10000,histtype='step',linewidth=1,linestyle='dashdot',color='k',range=[-1,1.01],density=True,cumulative=True,weights=a0all[sel]**(-1.)) line4, = ax.plot([0,1],[2,2],linewidth=1,linestyle='dashdot',color='k') sel = np.where((qall>8)&(qall<12)&(npassall[23,:]<6)&(a0all<2.5e4))[0] ax.hist(np.sin(2.*omegaall[sel]),bins=10000,histtype='step',linewidth=1,linestyle=(0,(3,5,1,5,1,5)),color='k',range=[-1,1.01],density=True,cumulative=True,weights=a0all[sel]**(-1.)) line5, = ax.plot([0,1],[2,2],linewidth=1,linestyle=(0,(3,5,1,5,1,5)),color='k') lineall = [line1, line2, line3, line4, line5] laball = ['Observed',r'N$_{\rm pass}$ = $\infty$',r'N$_{\rm pass}$ = 5',r'N$_{\rm pass}$ = 5, a$_0$ < 40,000 au',r'N$_{\rm pass}$ = 5, a$_0$ < 25,000 au'] ax.set_xlabel(r'$\sin$ 2$\omega_G$') ax.set_ylabel('Cumulative Fraction') ax.set_xlim(-1,1) ax.set_ylim(0,1) ax.axvspan(-1,0,alpha=0.2,color='k') ax.text(-.7,.9,'Increasing\nPerihelion',fontsize=18) ax.text(.3,.9,'Decreasing\nPerihelion',fontsize=18) ax.text(-.95,.92,'A',fontsize=16) ax.legend(lineall,laball,loc='center left') ax = fig.add_subplot(2,1,2) amaxvals = np.arange(2.5e4,5.1e4,1e3) testpvals = np.zeros(len(amaxvals)) for i,amaxval in enumerate(amaxvals): sel = np.where((npassall[23,:]<6)&(qall>8)&(qall<12)&(a0all8)&(qall<12)&(npassall[23,:]<2e5))[0] #ax.hist(np.sin(2.*omegaall[sel]),bins=10000,histtype='step',linewidth=1,linestyle='dotted',color='k',range=[-1,1.01],density=True,cumulative=True) #line2, = ax.plot([0,1],[2,2],linewidth=1,linestyle='dotted',color='k') #sel = np.where((qall>8)&(qall<12)&(npassall[23,:]<6))[0] #ax.hist(np.sin(2.*omegaall[sel]),bins=10000,histtype='step',linewidth=1,linestyle='dashed',color='k',range=[-1,1.01],density=True,cumulative=True) #line3, = ax.plot([0,1],[2,2],linewidth=1,linestyle='dashed',color='k') #sel = np.where((qall>8)&(qall<12)&(npassall[23,:]<6)&(a0all<4e4))[0] #ax.hist(np.sin(2.*omegaall[sel]),bins=10000,histtype='step',linewidth=1,linestyle='dashdot',color='k',range=[-1,1.01],density=True,cumulative=True,weights=a0all[sel]**(-1.)) #line4, = ax.plot([0,1],[2,2],linewidth=1,linestyle='dashdot',color='k') #sel = np.where((qall>8)&(qall<12)&(npassall[23,:]<6)&(a0all<2.5e4))[0] #ax.hist(np.sin(2.*omegaall[sel]),bins=10000,histtype='step',linewidth=1,linestyle=(0,(3,5,1,5,1,5)),color='k',range=[-1,1.01],density=True,cumulative=True,weights=a0all[sel]**(-1.)) #line5, = ax.plot([0,1],[2,2],linewidth=1,linestyle=(0,(3,5,1,5,1,5)),color='k') #lineall = [line1, line2, line3, line4, line5] #laball = ['Observed',r'N$_{\rm pass}$ = $\infty$',r'N$_{\rm pass}$ = 5',r'N$_{\rm pass}$ = 5, a$_0$ < 40,000 au',r'N$_{\rm pass}$ = 5, a$_0$ < 25,000 au'] ax.set_xlabel(r'$\sin$ 2$\omega_G$') ax.set_ylabel('Cumulative Fraction') ax.set_xlim(-1,1) ax.set_ylim(0,1) #ax.axvspan(-1,0,alpha=0.2,color='k') ax.text(-.7,.9,'Increasing\nPerihelion',fontsize=18) ax.text(.3,.9,'Decreasing\nPerihelion',fontsize=18) #ax.text(-.95,.92,'A',fontsize=16) #ax.legend(lineall,laball,loc='center left') fig.tight_layout() canvas.print_figure(figpath) figpath = 'OmegaMatch_talk2.pdf' fig = Figure(figsize=(8,6)) canvas = FigureCanvas(fig) ax = fig.add_subplot(1,1,1) ax.hist(np.sin(omegag[0:]*2.*math.pi/180.),bins=10000,histtype='step',linewidth=5,color='k',range=[-1,1.01],density=True,cumulative=True,alpha=0.5) line1, = ax.plot([0,1],[2,2],linewidth=5,color='k',alpha=0.5) sel = np.where((qall>8)&(qall<12)&(npassall[23,:]<2e5))[0] ax.hist(np.sin(2.*omegaall[sel]),bins=10000,histtype='step',linewidth=1,linestyle='dotted',color='k',range=[-1,1.01],density=True,cumulative=True) line2, = ax.plot([0,1],[2,2],linewidth=1,linestyle='dotted',color='k') #sel = np.where((qall>8)&(qall<12)&(npassall[23,:]<6))[0] #ax.hist(np.sin(2.*omegaall[sel]),bins=10000,histtype='step',linewidth=1,linestyle='dashed',color='k',range=[-1,1.01],density=True,cumulative=True) #line3, = ax.plot([0,1],[2,2],linewidth=1,linestyle='dashed',color='k') #sel = np.where((qall>8)&(qall<12)&(npassall[23,:]<6)&(a0all<4e4))[0] #ax.hist(np.sin(2.*omegaall[sel]),bins=10000,histtype='step',linewidth=1,linestyle='dashdot',color='k',range=[-1,1.01],density=True,cumulative=True,weights=a0all[sel]**(-1.)) #line4, = ax.plot([0,1],[2,2],linewidth=1,linestyle='dashdot',color='k') #sel = np.where((qall>8)&(qall<12)&(npassall[23,:]<6)&(a0all<2.5e4))[0] #ax.hist(np.sin(2.*omegaall[sel]),bins=10000,histtype='step',linewidth=1,linestyle=(0,(3,5,1,5,1,5)),color='k',range=[-1,1.01],density=True,cumulative=True,weights=a0all[sel]**(-1.)) #line5, = ax.plot([0,1],[2,2],linewidth=1,linestyle=(0,(3,5,1,5,1,5)),color='k') #lineall = [line1, line2, line3, line4, line5] #laball = ['Observed',r'N$_{\rm pass}$ = $\infty$',r'N$_{\rm pass}$ = 5',r'N$_{\rm pass}$ = 5, a$_0$ < 40,000 au',r'N$_{\rm pass}$ = 5, a$_0$ < 25,000 au'] lineall = [line1, line2] laball = ['Observed',r'N$_{\rm pass}$ = $\infty$'] ax.set_xlabel(r'$\sin$ 2$\omega_G$') ax.set_ylabel('Cumulative Fraction') ax.set_xlim(-1,1) ax.set_ylim(0,1) #ax.axvspan(-1,0,alpha=0.2,color='k') ax.text(-.7,.9,'Increasing\nPerihelion',fontsize=18) ax.text(.3,.9,'Decreasing\nPerihelion',fontsize=18) #ax.text(-.95,.92,'A',fontsize=16) ax.legend(lineall,laball,loc='center left') fig.tight_layout() canvas.print_figure(figpath) figpath = 'OmegaMatch_talk3.pdf' fig = Figure(figsize=(8,6)) canvas = FigureCanvas(fig) ax = fig.add_subplot(1,1,1) ax.hist(np.sin(omegag[0:]*2.*math.pi/180.),bins=10000,histtype='step',linewidth=5,color='k',range=[-1,1.01],density=True,cumulative=True,alpha=0.5) line1, = ax.plot([0,1],[2,2],linewidth=5,color='k',alpha=0.5) sel = np.where((qall>8)&(qall<12)&(npassall[23,:]<2e5))[0] ax.hist(np.sin(2.*omegaall[sel]),bins=10000,histtype='step',linewidth=1,linestyle='dotted',color='k',range=[-1,1.01],density=True,cumulative=True) line2, = ax.plot([0,1],[2,2],linewidth=1,linestyle='dotted',color='k') sel = np.where((qall>8)&(qall<12)&(npassall[23,:]<6))[0] ax.hist(np.sin(2.*omegaall[sel]),bins=10000,histtype='step',linewidth=1,linestyle='dashed',color='k',range=[-1,1.01],density=True,cumulative=True) line3, = ax.plot([0,1],[2,2],linewidth=1,linestyle='dashed',color='k') #sel = np.where((qall>8)&(qall<12)&(npassall[23,:]<6)&(a0all<4e4))[0] #ax.hist(np.sin(2.*omegaall[sel]),bins=10000,histtype='step',linewidth=1,linestyle='dashdot',color='k',range=[-1,1.01],density=True,cumulative=True,weights=a0all[sel]**(-1.)) #line4, = ax.plot([0,1],[2,2],linewidth=1,linestyle='dashdot',color='k') #sel = np.where((qall>8)&(qall<12)&(npassall[23,:]<6)&(a0all<2.5e4))[0] #ax.hist(np.sin(2.*omegaall[sel]),bins=10000,histtype='step',linewidth=1,linestyle=(0,(3,5,1,5,1,5)),color='k',range=[-1,1.01],density=True,cumulative=True,weights=a0all[sel]**(-1.)) #line5, = ax.plot([0,1],[2,2],linewidth=1,linestyle=(0,(3,5,1,5,1,5)),color='k') #lineall = [line1, line2, line3, line4, line5] #laball = ['Observed',r'N$_{\rm pass}$ = $\infty$',r'N$_{\rm pass}$ = 5',r'N$_{\rm pass}$ = 5, a$_0$ < 40,000 au',r'N$_{\rm pass}$ = 5, a$_0$ < 25,000 au'] lineall = [line1, line2, line3] laball = ['Observed',r'N$_{\rm pass}$ = $\infty$',r'N$_{\rm pass}$ = 5'] ax.set_xlabel(r'$\sin$ 2$\omega_G$') ax.set_ylabel('Cumulative Fraction') ax.set_xlim(-1,1) ax.set_ylim(0,1) #ax.axvspan(-1,0,alpha=0.2,color='k') ax.text(-.7,.9,'Increasing\nPerihelion',fontsize=18) ax.text(.3,.9,'Decreasing\nPerihelion',fontsize=18) #ax.text(-.95,.92,'A',fontsize=16) ax.legend(lineall,laball,loc='center left') fig.tight_layout() canvas.print_figure(figpath) figpath = 'OmegaMatch_talk4.pdf' fig = Figure(figsize=(8,6)) canvas = FigureCanvas(fig) ax = fig.add_subplot(1,1,1) ax.hist(np.sin(omegag[0:]*2.*math.pi/180.),bins=10000,histtype='step',linewidth=5,color='k',range=[-1,1.01],density=True,cumulative=True,alpha=0.5) line1, = ax.plot([0,1],[2,2],linewidth=5,color='k',alpha=0.5) sel = np.where((qall>8)&(qall<12)&(npassall[23,:]<2e5))[0] ax.hist(np.sin(2.*omegaall[sel]),bins=10000,histtype='step',linewidth=1,linestyle='dotted',color='k',range=[-1,1.01],density=True,cumulative=True) line2, = ax.plot([0,1],[2,2],linewidth=1,linestyle='dotted',color='k') sel = np.where((qall>8)&(qall<12)&(npassall[23,:]<6))[0] ax.hist(np.sin(2.*omegaall[sel]),bins=10000,histtype='step',linewidth=1,linestyle='dashed',color='k',range=[-1,1.01],density=True,cumulative=True) line3, = ax.plot([0,1],[2,2],linewidth=1,linestyle='dashed',color='k') sel = np.where((qall>8)&(qall<12)&(npassall[23,:]<6)&(a0all<4e4))[0] ax.hist(np.sin(2.*omegaall[sel]),bins=10000,histtype='step',linewidth=1,linestyle='dashdot',color='k',range=[-1,1.01],density=True,cumulative=True,weights=a0all[sel]**(-1.)) line4, = ax.plot([0,1],[2,2],linewidth=1,linestyle='dashdot',color='k') #sel = np.where((qall>8)&(qall<12)&(npassall[23,:]<6)&(a0all<2.5e4))[0] #ax.hist(np.sin(2.*omegaall[sel]),bins=10000,histtype='step',linewidth=1,linestyle=(0,(3,5,1,5,1,5)),color='k',range=[-1,1.01],density=True,cumulative=True,weights=a0all[sel]**(-1.)) #line5, = ax.plot([0,1],[2,2],linewidth=1,linestyle=(0,(3,5,1,5,1,5)),color='k') #lineall = [line1, line2, line3, line4, line5] #laball = ['Observed',r'N$_{\rm pass}$ = $\infty$',r'N$_{\rm pass}$ = 5',r'N$_{\rm pass}$ = 5, a$_0$ < 40,000 au',r'N$_{\rm pass}$ = 5, a$_0$ < 25,000 au'] lineall = [line1, line2, line3, line4] laball = ['Observed',r'N$_{\rm pass}$ = $\infty$',r'N$_{\rm pass}$ = 5',r'N$_{\rm pass}$ = 5, a$_0$ < 40,000 au'] ax.set_xlabel(r'$\sin$ 2$\omega_G$') ax.set_ylabel('Cumulative Fraction') ax.set_xlim(-1,1) ax.set_ylim(0,1) #ax.axvspan(-1,0,alpha=0.2,color='k') ax.text(-.7,.9,'Increasing\nPerihelion',fontsize=18) ax.text(.3,.9,'Decreasing\nPerihelion',fontsize=18) #ax.text(-.95,.92,'A',fontsize=16) ax.legend(lineall,laball,loc='center left') fig.tight_layout() canvas.print_figure(figpath) figpath = 'OmegaMatch_talk5.pdf' fig = Figure(figsize=(8,6)) canvas = FigureCanvas(fig) ax = fig.add_subplot(1,1,1) ax.hist(np.sin(omegag[0:]*2.*math.pi/180.),bins=10000,histtype='step',linewidth=5,color='k',range=[-1,1.01],density=True,cumulative=True,alpha=0.5) line1, = ax.plot([0,1],[2,2],linewidth=5,color='k',alpha=0.5) sel = np.where((qall>8)&(qall<12)&(npassall[23,:]<2e5))[0] ax.hist(np.sin(2.*omegaall[sel]),bins=10000,histtype='step',linewidth=1,linestyle='dotted',color='k',range=[-1,1.01],density=True,cumulative=True) line2, = ax.plot([0,1],[2,2],linewidth=1,linestyle='dotted',color='k') sel = np.where((qall>8)&(qall<12)&(npassall[23,:]<6))[0] ax.hist(np.sin(2.*omegaall[sel]),bins=10000,histtype='step',linewidth=1,linestyle='dashed',color='k',range=[-1,1.01],density=True,cumulative=True) line3, = ax.plot([0,1],[2,2],linewidth=1,linestyle='dashed',color='k') sel = np.where((qall>8)&(qall<12)&(npassall[23,:]<6)&(a0all<4e4))[0] ax.hist(np.sin(2.*omegaall[sel]),bins=10000,histtype='step',linewidth=1,linestyle='dashdot',color='k',range=[-1,1.01],density=True,cumulative=True,weights=a0all[sel]**(-1.)) line4, = ax.plot([0,1],[2,2],linewidth=1,linestyle='dashdot',color='k') sel = np.where((qall>8)&(qall<12)&(npassall[23,:]<6)&(a0all<2.5e4))[0] ax.hist(np.sin(2.*omegaall[sel]),bins=10000,histtype='step',linewidth=1,linestyle=(0,(3,5,1,5,1,5)),color='k',range=[-1,1.01],density=True,cumulative=True,weights=a0all[sel]**(-1.)) line5, = ax.plot([0,1],[2,2],linewidth=1,linestyle=(0,(3,5,1,5,1,5)),color='k') lineall = [line1, line2, line3, line4, line5] laball = ['Observed',r'N$_{\rm pass}$ = $\infty$',r'N$_{\rm pass}$ = 5',r'N$_{\rm pass}$ = 5, a$_0$ < 40,000 au',r'N$_{\rm pass}$ = 5, a$_0$ < 25,000 au'] #lineall = [line1, line2, line3, line4] #laball = ['Observed',r'N$_{\rm pass}$ = $\infty$',r'N$_{\rm pass}$ = 5','N$_{\rm pass}$ = 5, a$_0$ < 40,000 au'] ax.set_xlabel(r'$\sin$ 2$\omega_G$') ax.set_ylabel('Cumulative Fraction') ax.set_xlim(-1,1) ax.set_ylim(0,1) #ax.axvspan(-1,0,alpha=0.2,color='k') ax.text(-.7,.9,'Increasing\nPerihelion',fontsize=18) ax.text(.3,.9,'Decreasing\nPerihelion',fontsize=18) #ax.text(-.95,.92,'A',fontsize=16) ax.legend(lineall,laball,loc='center left') fig.tight_layout() canvas.print_figure(figpath) figpath = 'OmegaMatch_talk6.pdf' fig = Figure(figsize=(8,6)) canvas = FigureCanvas(fig) ax = fig.add_subplot(1,1,1) amaxvals = np.arange(2.5e4,5.1e4,1e3) testpvals = np.zeros(len(amaxvals)) for i,amaxval in enumerate(amaxvals): sel = np.where((npassall[23,:]<6)&(qall>8)&(qall<12)&(a0all8))[0] ax.hist(qprevall[sel],bins=100000,range=[0,10000],histtype='step',cumulative=True,density=True,linewidth=1,color='k',linestyle='dotted') line2, = ax.plot([0,1],[2,2],linewidth=1,color='k',linestyle='dotted') sel = np.where((npassall[23,:]<6)&(qall<12)&(qall>8))[0] ax.hist(qprevall[sel],bins=100000,range=[0,10000],histtype='step',cumulative=True,density=True,linewidth=1,color='k',linestyle='dashed') line3, = ax.plot([0,1],[2,2],linewidth=1,color='k',linestyle='dashed') #sel = np.where((npassall[23,:]<6)&(qall<12)&(qall>8)&(a0all<4e4))[0] #ax.hist(qprevall[sel],bins=100000,range=[0,10000],histtype='step',cumulative=True,density=True,linewidth=1,color='k',linestyle='solid') #line3, = ax.plot([0,1],[2,2],linewidth=1,color='k',linestyle='dashed') lineall = [line1,line2,line3] laball = ['Observed',r'N$_{\rm Pass}$ = $\infty$',r'N$_{\rm Pass}$ = 5'] ax.set_xlim(0,32) ax.set_ylim(0,1) ax.set_xlabel('Previous Perihelion Passage (au)') ax.set_ylabel('Cumulative Fraction') ax.legend(lineall,laball,loc='upper left') ax.text(5.,.03,'J',fontsize=18) ax.text(9.3,.03,'S',fontsize=18) ax.text(19.1,.03,'U',fontsize=18) ax.text(29.9,.03,'N',fontsize=18) fig.tight_layout() canvas.print_figure(figpath) figpath = 'qprevs_talk1.pdf' fig = Figure(figsize=(8,6)) canvas = FigureCanvas(fig) ax = fig.add_subplot(1,1,1) ax.hist(qprevreal,bins=100000,range=[0,10000],histtype='step',cumulative=True,density=True,linewidth=3,color='k',alpha=0.5,linestyle='solid') line1, = ax.plot([0,1],[2,2],linewidth=3,alpha=0.5,color='k',linestyle='solid') #sel = np.where((npassall[23,:]<1e5)&(qall<12)&(qall>8))[0] #ax.hist(qprevall[sel],bins=100000,range=[0,10000],histtype='step',cumulative=True,density=True,linewidth=1,color='r',linestyle='dotted') #line2, = ax.plot([0,1],[2,2],linewidth=1,color='r',linestyle='dotted') #sel = np.where((npassall[23,:]<6)&(qall<12)&(qall>8))[0] #ax.hist(qprevall[sel],bins=100000,range=[0,10000],histtype='step',cumulative=True,density=True,linewidth=1,color='b',linestyle='dashed') #line3, = ax.plot([0,1],[2,2],linewidth=1,color='b',linestyle='dashed') #sel = np.where((npassall[23,:]<6)&(qall<12)&(qall>8)&(a0all<4e4))[0] #ax.hist(qprevall[sel],bins=100000,range=[0,10000],histtype='step',cumulative=True,density=True,linewidth=1,color='k',linestyle='solid') #line3, = ax.plot([0,1],[2,2],linewidth=1,color='k',linestyle='dashed') lineall = [line1,line2,line3] laball = ['Observed',r'N$_{\rm Pass}$ = $\infty$',r'N$_{\rm Pass}$ = 5'] ax.set_xlim(0,32) ax.set_ylim(0,1) ax.set_xlabel('Previous Perihelion Passage (au)') ax.set_ylabel('Cumulative Fraction') #ax.legend(lineall,laball,loc='upper left') ax.text(5.,.03,'J',fontsize=18) ax.text(9.3,.03,'S',fontsize=18) ax.text(19.1,.03,'U',fontsize=18) ax.text(29.9,.03,'N',fontsize=18) fig.tight_layout() canvas.print_figure(figpath) figpath = 'qprevs_talk2.pdf' fig = Figure(figsize=(8,6)) canvas = FigureCanvas(fig) ax = fig.add_subplot(1,1,1) ax.hist(qprevreal,bins=100000,range=[0,10000],histtype='step',cumulative=True,density=True,linewidth=3,color='k',alpha=0.5,linestyle='solid') line1, = ax.plot([0,1],[2,2],linewidth=3,alpha=0.5,color='k',linestyle='solid') #sel = np.where((npassall[23,:]<1e5)&(qall<12)&(qall>8))[0] #ax.hist(qprevall[sel],bins=100000,range=[0,10000],histtype='step',cumulative=True,density=True,linewidth=1,color='r',linestyle='dotted') #line2, = ax.plot([0,1],[2,2],linewidth=1,color='r',linestyle='dotted') sel = np.where((npassall[23,:]<6)&(qall<12)&(qall>8))[0] ax.hist(qprevall[sel],bins=100000,range=[0,10000],histtype='step',cumulative=True,density=True,linewidth=1,color='b',linestyle='dashed') line3, = ax.plot([0,1],[2,2],linewidth=1,color='b',linestyle='dashed') #sel = np.where((npassall[23,:]<6)&(qall<12)&(qall>8)&(a0all<4e4))[0] #ax.hist(qprevall[sel],bins=100000,range=[0,10000],histtype='step',cumulative=True,density=True,linewidth=1,color='k',linestyle='solid') #line3, = ax.plot([0,1],[2,2],linewidth=1,color='k',linestyle='dashed') lineall = [line1,line3] laball = ['Observed',r'N$_{\rm Pass}$ = 5'] ax.set_xlim(0,32) ax.set_ylim(0,1) ax.set_xlabel('Previous Perihelion Passage (au)') ax.set_ylabel('Cumulative Fraction') ax.legend(lineall,laball,loc='upper left') ax.text(5.,.03,'J',fontsize=18) ax.text(9.3,.03,'S',fontsize=18) ax.text(19.1,.03,'U',fontsize=18) ax.text(29.9,.03,'N',fontsize=18) fig.tight_layout() canvas.print_figure(figpath) figpath = 'contourcombine_color.pdf' fig = Figure(figsize=(16,6)) canvas = FigureCanvas(fig) ax = fig.add_subplot(1,2,1) #ax.contourf(qcritvals,npassvals,apvals,[.0455,1],colors='b',alpha=0.2) #ax.contourf(qcritvals,npassvals,omegapvals,[.0455,1],colors='k',alpha=0.2) #ax.contourf(qcritvals,npassvals,qprevpvals,[.0455,1],colors='g',alpha=0.2) ax.contourf(qcritvals,npassvals,apvals,[.0455,1],colors='#377eb8',alpha=0.2) ax.contourf(qcritvals,npassvals,omegapvals,[.0455,1],colors='#999999',alpha=0.2) ax.contourf(qcritvals,npassvals,qprevpvals,[.0455,1],colors='#f781bf',alpha=0.2) #c1 = ax.contour(qcritvals,npassvals,apvals,[.0455,1],colors='b') c1 = ax.contour(qcritvals,npassvals,apvals,[.0455,1],colors='#377eb8') fmt1 = {} strs = ['a','a'] for l, s in zip(c1.levels, strs): fmt1[l] = s #c2 = ax.contour(qcritvals,npassvals,omegapvals,[.0455,1],colors='k') c2 = ax.contour(qcritvals,npassvals,omegapvals,[.0455,1],colors='#999999') fmt2 = {} strs = [r'sin 2$\omega_g$',r'sin 2$\omega_g$'] for l, s in zip(c2.levels, strs): fmt2[l] = s #c3 = ax.contour(qcritvals,npassvals,qprevpvals,[.0455,1],colors='g') c3 = ax.contour(qcritvals,npassvals,qprevpvals,[.0455,1],colors='#f781bf') fmt3 = {} strs = [r'q$_{\rm prev}$',r'q$_{\rm prev}$'] for l, s in zip(c3.levels, strs): fmt3[l] = s ax.clabel(c1,c1.levels,inline=True,fontsize=12,fmt=fmt1,inline_spacing=2) #ax.clabel(c1,c1.levels,inline=true,fontsize=10,fmt=fmt1) ax.clabel(c2,c2.levels,inline=True,fontsize=10,fmt=fmt2,inline_spacing=1) #ax.clabel(c2,c2.levels,inline=true,fontsize=10,fmt=fmt2) ax.clabel(c3,c3.levels,inline=True,fontsize=11,fmt=fmt3,inline_spacing=1) #ax.clabel(c3,c3.levels,inline=true,fontsize=10,fmt=fmt3) #add near-Earth LPC percentiles to the plot near = np.where((qnear>0)&(qnear<4)&(anear>1e4))[0] npass50per = np.percentile(npassnear[:,near],50,axis=1) npass75per = np.percentile(npassnear[:,near],80,axis=1) npass90per = np.percentile(npassnear[:,near],90,axis=1) #ax.plot(qcritvals,npass50per,color='k',linewidth=1,linestyle='dotted') #ax.plot(qcritvals,npass75per,color='k',linewidth=1,linestyle='dotted') ax.plot(qcritvals,npass75per-1,color='k',linewidth=2,linestyle='dashed') ax.set_ylim(1,14) ax.set_xlim(10,20) ax.set_xlabel(r'q$_{\rm Fade}$ (au)',fontsize=16) ax.set_ylabel(r'N$_{\rm Pass}$',fontsize=16) ax.set_title(r'20,000 < a$_{\rm Initial}$ < 50,000 au') ax.text(10.2, 13.3, 'A', fontsize=18) ax.text(13.5,2.1,'Near-Earth LPCs 80th Percentile',rotation=0) ax = fig.add_subplot(1,2,2) ax.contourf(qcritvals,npassvals,apvals_acap,[.0455,1],colors='#377eb8',alpha=0.2) ax.contourf(qcritvals,npassvals,omegapvals_acap,[.0455,1],colors='#999999',alpha=0.2) ax.contourf(qcritvals,npassvals,qprevpvals_acap,[.0455,1],colors='#f781bf',alpha=0.2) c1 = ax.contour(qcritvals,npassvals,apvals_acap,[.0455,1],colors='#377eb8') fmt1 = {} strs = ['a','a'] for l, s in zip(c1.levels, strs): fmt1[l] = s c2 = ax.contour(qcritvals,npassvals,omegapvals_acap,[.0455,1],colors='#999999') fmt2 = {} strs = [r'sin 2$\omega_g$',r'sin 2$\omega_g$'] for l, s in zip(c2.levels, strs): fmt2[l] = s c3 = ax.contour(qcritvals,npassvals,qprevpvals_acap,[.0455,1],colors='#f781bf') fmt3 = {} strs = [r'q$_{\rm prev}$',r'q$_{\rm prev}$'] for l, s in zip(c3.levels, strs): fmt3[l] = s ax.clabel(c1,c1.levels,inline=True,fontsize=10,fmt=fmt1,inline_spacing=2) #ax.clabel(c1,c1.levels,inline=true,fontsize=10,fmt=fmt1) ax.clabel(c2,c2.levels,inline=True,fontsize=12,fmt=fmt2,inline_spacing=0) #ax.clabel(c2,c2.levels,inline=true,fontsize=10,fmt=fmt2) ax.clabel(c3,c3.levels,inline=True,fontsize=11,fmt=fmt3,inline_spacing=1) #ax.clabel(c3,c3.levels,inline=true,fontsize=10,fmt=fmt3) #add near-Earth LPC percentiles to the plot near = np.where((qnear>0)&(qnear<4)&(anear>1e4)&(a0near<2.5e4))[0] npass50per = np.percentile(npassnear[:,near],50,axis=1) npass75per = np.percentile(npassnear[:,near],80,axis=1) npass90per = np.percentile(npassnear[:,near],90,axis=1) #ax.plot(qcritvals,npass50per,color='k',linewidth=1,linestyle='dotted') #ax.plot(qcritvals,npass75per,color='k',linewidth=1,linestyle='dotted') ax.plot(qcritvals,npass75per-1,color='k',linewidth=2,linestyle='dashed') ax.set_ylim(1,14) ax.set_xlim(10,20) ax.set_xlabel(r'q$_{\rm Fade}$ (au)',fontsize=16) ax.set_ylabel(r'N$_{\rm Pass}$',fontsize=16) ax.set_title(r'20,000 < a$_{\rm Initial}$ < 25,000 au') ax.text(10.2, 13.3, 'B', fontsize=18) ax.text(14.8,4.0,'Near-Earth LPCs 80th Percentile',rotation=12) fig.tight_layout() canvas.print_figure(figpath) figpath = 'contourcombine_color_talk1.pdf' fig = Figure(figsize=(8,6)) canvas = FigureCanvas(fig) ax = fig.add_subplot(1,1,1) #ax.contourf(qcritvals,npassvals,apvals,[.0455,1],colors='b',alpha=0.2) #ax.contourf(qcritvals,npassvals,omegapvals,[.0455,1],colors='k',alpha=0.2) #ax.contourf(qcritvals,npassvals,qprevpvals,[.0455,1],colors='g',alpha=0.2) ax.contourf(qcritvals,npassvals,apvals,[.0455,1],colors='#377eb8',alpha=0.2) ax.contourf(qcritvals,npassvals,omegapvals,[.0455,1],colors='#999999',alpha=0.2) ax.contourf(qcritvals,npassvals,qprevpvals,[.0455,1],colors='#f781bf',alpha=0.2) #c1 = ax.contour(qcritvals,npassvals,apvals,[.0455,1],colors='b') c1 = ax.contour(qcritvals,npassvals,apvals,[.0455,1],colors='#377eb8') fmt1 = {} strs = ['a','a'] for l, s in zip(c1.levels, strs): fmt1[l] = s #c2 = ax.contour(qcritvals,npassvals,omegapvals,[.0455,1],colors='k') c2 = ax.contour(qcritvals,npassvals,omegapvals,[.0455,1],colors='#999999') fmt2 = {} strs = [r'sin 2$\omega_g$',r'sin 2$\omega_g$'] for l, s in zip(c2.levels, strs): fmt2[l] = s #c3 = ax.contour(qcritvals,npassvals,qprevpvals,[.0455,1],colors='g') c3 = ax.contour(qcritvals,npassvals,qprevpvals,[.0455,1],colors='#f781bf') fmt3 = {} strs = [r'q$_{\rm prev}$',r'q$_{\rm prev}$'] for l, s in zip(c3.levels, strs): fmt3[l] = s ax.clabel(c1,c1.levels,inline=True,fontsize=12,fmt=fmt1,inline_spacing=2) #ax.clabel(c1,c1.levels,inline=true,fontsize=10,fmt=fmt1) ax.clabel(c2,c2.levels,inline=True,fontsize=10,fmt=fmt2,inline_spacing=1) #ax.clabel(c2,c2.levels,inline=true,fontsize=10,fmt=fmt2) ax.clabel(c3,c3.levels,inline=True,fontsize=11,fmt=fmt3,inline_spacing=1) #ax.clabel(c3,c3.levels,inline=true,fontsize=10,fmt=fmt3) #add near-Earth LPC percentiles to the plot near = np.where((qnear>0)&(qnear<4)&(anear>1e4))[0] npass50per = np.percentile(npassnear[:,near],50,axis=1) npass75per = np.percentile(npassnear[:,near],80,axis=1) npass90per = np.percentile(npassnear[:,near],90,axis=1) #ax.plot(qcritvals,npass50per,color='k',linewidth=1,linestyle='dotted') #ax.plot(qcritvals,npass75per,color='k',linewidth=1,linestyle='dotted') #ax.plot(qcritvals,npass75per-1,color='k',linewidth=2,linestyle='dashed') ax.set_ylim(1,14) ax.set_xlim(10,20) ax.set_xlabel(r'q$_{\rm Fade}$ (au)',fontsize=16) ax.set_ylabel(r'N$_{\rm Pass}$',fontsize=16) #ax.set_title(r'20,000 < a$_{\rm Initial}$ < 50,000 au') #ax.text(10.2, 13.3, 'A', fontsize=18) #ax.text(13.5,2.1,'Near-Earth LPCs 80th Percentile',rotation=0) fig.tight_layout() canvas.print_figure(figpath) figpath = 'contourcombine_color_talk2.pdf' fig = Figure(figsize=(8,6)) canvas = FigureCanvas(fig) ax = fig.add_subplot(1,1,1) #ax.contourf(qcritvals,npassvals,apvals,[.0455,1],colors='b',alpha=0.2) #ax.contourf(qcritvals,npassvals,omegapvals,[.0455,1],colors='k',alpha=0.2) #ax.contourf(qcritvals,npassvals,qprevpvals,[.0455,1],colors='g',alpha=0.2) ax.contourf(qcritvals,npassvals,apvals,[.0455,1],colors='#377eb8',alpha=0.2) ax.contourf(qcritvals,npassvals,omegapvals,[.0455,1],colors='#999999',alpha=0.2) ax.contourf(qcritvals,npassvals,qprevpvals,[.0455,1],colors='#f781bf',alpha=0.2) #c1 = ax.contour(qcritvals,npassvals,apvals,[.0455,1],colors='b') c1 = ax.contour(qcritvals,npassvals,apvals,[.0455,1],colors='#377eb8') fmt1 = {} strs = ['a','a'] for l, s in zip(c1.levels, strs): fmt1[l] = s #c2 = ax.contour(qcritvals,npassvals,omegapvals,[.0455,1],colors='k') c2 = ax.contour(qcritvals,npassvals,omegapvals,[.0455,1],colors='#999999') fmt2 = {} strs = [r'sin 2$\omega_g$',r'sin 2$\omega_g$'] for l, s in zip(c2.levels, strs): fmt2[l] = s #c3 = ax.contour(qcritvals,npassvals,qprevpvals,[.0455,1],colors='g') c3 = ax.contour(qcritvals,npassvals,qprevpvals,[.0455,1],colors='#f781bf') fmt3 = {} strs = [r'q$_{\rm prev}$',r'q$_{\rm prev}$'] for l, s in zip(c3.levels, strs): fmt3[l] = s ax.clabel(c1,c1.levels,inline=True,fontsize=12,fmt=fmt1,inline_spacing=2) #ax.clabel(c1,c1.levels,inline=true,fontsize=10,fmt=fmt1) ax.clabel(c2,c2.levels,inline=True,fontsize=10,fmt=fmt2,inline_spacing=1) #ax.clabel(c2,c2.levels,inline=true,fontsize=10,fmt=fmt2) ax.clabel(c3,c3.levels,inline=True,fontsize=11,fmt=fmt3,inline_spacing=1) #ax.clabel(c3,c3.levels,inline=true,fontsize=10,fmt=fmt3) #add near-Earth LPC percentiles to the plot near = np.where((qnear>0)&(qnear<4)&(anear>1e4))[0] npass50per = np.percentile(npassnear[:,near],50,axis=1) npass75per = np.percentile(npassnear[:,near],80,axis=1) npass90per = np.percentile(npassnear[:,near],90,axis=1) #ax.plot(qcritvals,npass50per,color='k',linewidth=1,linestyle='dotted') #ax.plot(qcritvals,npass75per,color='k',linewidth=1,linestyle='dotted') ax.plot(qcritvals,npass75per-1,color='k',linewidth=2,linestyle='dashed') ax.set_ylim(1,14) ax.set_xlim(10,20) ax.set_xlabel(r'q$_{\rm Fade}$ (au)',fontsize=16) ax.set_ylabel(r'N$_{\rm Pass}$',fontsize=16) #ax.set_title(r'20,000 < a$_{\rm Initial}$ < 50,000 au') #ax.text(10.2, 13.3, 'A', fontsize=18) ax.text(13.5,2.1,'Near-Earth LPCs 80th Percentile',rotation=0) fig.tight_layout() canvas.print_figure(figpath) cmd = 'say \"done now\"' os.system(cmd)