from const import * import numpy as np import matplotlib import matplotlib.pyplot as plt import sys, time ##################################################### # I'm defining functions here for rotation matrices # ##################################################### def r1(a): '''This returns a rotation by angle of a radians CCW around the x axis ''' return np.array([[1,0,0],[0, np.cos(a), -np.sin(a)],[0, np.sin(a), np.cos(a)]]) def r2(a): '''This returns a rotation by angle of a radians CCW around the y axis ''' return np.array([[np.cos(a),0,-np.sin(a)],[0, 1, 0],[np.sin(a), 0, np.cos(a)]]) def r3(a): '''This returns a rotation by angle of a radians CCW around the z axis ''' return np.array([[np.cos(a),-np.sin(a),0],[np.sin(a), np.cos(a), 0],[0, 0, 1]]) def norm(a): '''This returns the unit vector in the direction of a''' mag = np.linalg.norm(a) if mag == 0: print "Why the heck is |a| = 0?" return a return a / mag ######################################################################### ######################################################################### ############################ Pulsar vectors ############################ ######################################################################### Nhat = np.array([0,0,1]) # This is the northward vector # Unit vector from Earth to source rHat = np.array([np.cos(dlt) * np.cos(alp), np.cos(dlt) * np.sin(alp), np.sin(dlt)]) E = norm(np.cross(rHat, Nhat)) # Local E NLoc = norm(np.cross(E, rHat)) # Local North wx = norm(np.sin(psi) * NLoc + np.cos(psi) * E) # Rotating by psi wy = norm(np.cos(psi) * NLoc - np.sin(psi) * E) # Rotating by psi wz = -norm(rHat) # Z points toward Earth ######################################################################### class Curve(object): def __init__(self, h0, phi0, totalTime): self.strain = h0 self.angle = phi0 self.time = totalTime def lng(self, t): '''Returns the longitude of the detector in time; phi0 is angle at t = 0 relative that astronomical thing''' return 2 * np.pi * t / Td + self.angle def d_detector(self, t): '''This calculates the 3 3-vectors of the detector at time t''' longi= self.lng(t) dz = np.array([np.cos(lat) * np.cos(longi), np.cos(lat) * np.sin(longi), np.sin(lat)]) E = norm(np.cross(Nhat, dz)) # E is perp. plane of Nhat and dz NLoc = norm(np.cross(dz, E)) dx = norm(-np.sin(beta) * NLoc + np.cos(beta) * E) dy = norm(np.cos(beta) * NLoc + np.sin(beta) * E) dz *= rE return (dx, dy, dz) def APl(self, t): (dx, dy, dz) = self.d_detector(t) return (wx.dot(dx) ** 2 - wx.dot(dy) ** 2 - wy.dot(dx) ** 2 + wy.dot(dy) ** 2) / 2 def ACr(self, t): (dx, dy, dz) = self.d_detector(t) return (wx.dot(dx)) * (wy.dot(dx)) - (wx.dot(dy)) * (wy.dot(dy)) def makeAPs(self): tIn = np.arange(0, self.time, sampR) plus = np.zeros(len(tIn)) crss = np.zeros(len(tIn)) for i, e in enumerate(tIn): (dx, dy, dz) = self.d_detector(e) plus[i] = (wx.dot(dx) ** 2 - wx.dot(dy) ** 2 - wy.dot(dx) ** 2 + wy.dot(dy) ** 2) / 2 crss[i] = (wx.dot(dx)) * (wy.dot(dx)) - (wx.dot(dy)) * (wy.dot(dy)) return self.strain * (plus, crss) def generateCurves(self): (plus, crss) = self.makeAPs() # NOTE: These already include the strain return (1 + np.cos(iota) ** 2) / 4 * plus - 1j/2 * np.cos(iota) * crss class deltaTheta(object): def __init__(self, beta, tIn): self.beta = beta self.tIn = tIn def r_orbit(self, t): '''Gives rectangular coordinates of planet in orbit''' ang = 2 * np.pi * t.T / TYr # Angle over time; t is since a 0 angle pos = rOrb * np.array([[np.cos(ang)], [np.sin(ang)], [np.zeros(len(ang))]]) return pos.T.dot(r1(tilt)) # Because the Earth orbits at an angle relative N def r_revol(self, t): '''Gives rectangular coordinates of observatory relative center of Earth''' lngi = 2 * np.pi * t / Td # Angle over time; t is since a 0 in the angle return rE * np.array([[np.cos(lat) * np.cos(lngi)], [np.sin(lat) * np.cos(lngi)], [np.sin(lngi)]]).T def r(self, t): return self.r_revol(t) #+ self.r_orbit(t) def romerDelay(self, tIn): return self.r(tIn).dot(rHat) / c def expAngle0(self): ''' Old method to calculate expAngle; keeping in case''' l = len(self.tIn) r = np.zeros((l, 3)) for ind, t in enumerate(self.tIn): r[ind][0:3] = self.r(t).T k = r.dot(rHat) #rHat is like n hat dTheta_dT = [np.zeros(l), nu * np.ones(l) + nuD * self.tIn + nuDD / 2 * (self.tIn ** 2), nuD * np.ones(l) + nuDD * self.tIn, nuDD * np.ones(l)] factorial= [1, 1, 2, 6] rnge = np.arange(1, 4) delta = 0 for i in rnge: delta += (c ** -i) / (factorial[i]) * (self.beta ** -i -1) * (k ** i) * dTheta_dT[i] return np.exp(1j * delta) def expAngle(self, rD): '''Approximates nu dot = 0''' d = (self.beta ** -1 - 1) * nu * rD return np.exp(1j * d) if __name__ == "__main__": finT= Td * 200 tIn = np.arange(0, finT, sampR) # Base curve obj = Curve(1, 0, finT) gra = obj.generateCurves() # For creating the delay obj = deltaTheta(1.0001, tIn) roD = obj.romerDelay(tIn) add = obj.expAngle(roD) add = np.squeeze(add) gra2 = gra * add #plt.plot(gra.real) #plt.plot(gra.imag) plt.plot(gra2.real) plt.plot(gra2.imag) plt.show() FT = np.fft.fft(gra2) lamb2 = np.abs(FT) * np.sign(FT).real plt.plot(tIn / (sampR * finT), lamb2) plt.show() ''' # Times the generation of deltaTheta numTrials = 1000 timeLog = np.zeros(numTrials) for i in range(numTrials): t0 = time.clock() a = deltaTheta(float(sys.argv[1]), np.arange(0,Td,60)) b = a.expAngle() t1 = time.clock() timeLog[i] = t1 - t0 plt.figure() plt.hist(timeLog, numTrials, color="k", histtype="step") plt.title("Time To calculate") plt.show() #Currently derives distribution for time required to generate a curve. numTrials = 1000 timeLog = np.zeros(numTrials) for i in range(numTrials): t0 = time.clock() obj = Curve(float(sys.argv[1]), float(sys.argv[2]), float(sys.argv[3])) gra = obj.makeCurve() t1 = time.clock() timeLog[i] = t1 - t0 plt.figure() plt.hist(timeLog, numTrials / 5, color="k", histtype="step") plt.title("Time To calculate") plt.show() '''