Commit b73683f9 authored by Karsten Wiesner's avatar Karsten Wiesner 💬

Merge branch '10-add-long-arm-signal-calculation' into 'master'

Resolve "Add long arm signal calculation"

Closes #10

See merge request !5
parents 0c515785 9b8e1858
From 5ea4cc063d0522e1d88480a740e662a3f7516ff8 Mon Sep 17 00:00:00 2001
From: Martin Staab <martin.staab@aei.mpg.de>
Date: Fri, 7 Jun 2019 14:38:36 +0200
Subject: [PATCH 1/2] solved compatibility issues
---
fastsource/fastbinary/FastBinary.i | 21 +++++++++++----------
1 file changed, 11 insertions(+), 10 deletions(-)
diff --git a/fastsource/fastbinary/FastBinary.i b/fastsource/fastbinary/FastBinary.i
index b4d2349..7b2d287 100644
--- a/fastsource/fastbinary/FastBinary.i
+++ b/fastsource/fastbinary/FastBinary.i
@@ -28,10 +28,11 @@
%pythoncode %{
import sys, math, time
-if 'lisaxml2' in sys.modules:
- import lisaxml2 as lisaxml
-else:
- import lisaxml
+#if 'lisaxml2' in sys.modules:
+# import lisaxml2 as lisaxml
+#else:
+# import lisaxml
+import lisaxml2 as lisaxml
import numpy
import mc3lisa, FrequencyArray
@@ -69,7 +70,7 @@ class FastGalacticBinary(lisaxml.Source):
setattr(self,par[0],par[2])
def fourier(self,simulator='synthlisa',table=None,T=6.2914560e7,dt=15,algorithm='mldc',oversample=1,kmin=0,length=None,mpi=False,status=True):
- if table == None:
+ if table is None:
return self.onefourier(simulator=simulator,T=T,dt=dt,algorithm=algorithm,oversample=oversample)
elif mpi == False:
if length == None:
@@ -136,7 +137,7 @@ class FastGalacticBinary(lisaxml.Source):
MPI.COMM_WORLD.Recv((rbuf,MPI.DOUBLE_COMPLEX),s+1,i)
buf[i][:] += rbuf[:]
- # something strange going on with mpi4py's Reduce. Will use simple communication.
+ # something strange going on with mpi4pys Reduce. Will use simple communication.
# for i in range(3):
# MPI.COMM_WORLD.Reduce(None,(buf[i],MPI.DOUBLE_COMPLEX),MPI.SUM,0)
@@ -215,11 +216,11 @@ class FastGalacticBinary(lisaxml.Source):
Acut = A * math.sqrt(T/Sm)
# this is [2^(1+[log_2 SNR])]; remember that the response of the sinc is bound by 1/(delta bin #)
- # therefore we're using a number of bins proportional to SNR, rounded to the next power of two
+ # therefore we are using a number of bins proportional to SNR, rounded to the next power of two
# according to my tests, this is consistent with accepting +-0.4 in the detection SNR
M = int(math.pow(2.0,1 + int(math.log(Acut)/math.log(2.0))));
- # corrected per Neil's 2009-07-28 changes to Fast_Response3.c
+ # corrected per Neils 2009-07-28 changes to Fast_Response3.c
# use the larger of the two buffer sizes, but not above 8192
# -- further changed for new chirping logic (r1164)
if algorithm == 'mldc':
@@ -253,7 +254,7 @@ class FastGalacticBinary(lisaxml.Source):
return M,N
def onefourier(self,simulator='synthlisa',vector=None,buffer=None,T=6.2914560e7,dt=15,algorithm='mldc',oversample=1):
- if vector != None:
+ if not vector is None:
self.Frequency,self.FrequencyDerivative,self.Amplitude = vector[0],vector[1],vector[4]
M, N = self.buffersize(T,self.Frequency,self.FrequencyDerivative,self.Amplitude,algorithm,oversample)
@@ -272,7 +273,7 @@ class FastGalacticBinary(lisaxml.Source):
method = 0 if algorithm == 'legacy' else 1
# TO DO: pass a vector of parameters directly
- if vector != None:
+ if not vector is None:
fastbin.Response(vector[0],vector[1],0.5*math.pi - vector[2],vector[3],
vector[4],vector[5],vector[6],vector[7],
fastbin.XLS,fastbin.XSL,fastbin.YLS,fastbin.YSL,fastbin.ZLS,fastbin.ZSL,
--
2.21.0
From 22cb6723cf6df46b658b4ba1f926075456ff8a37 Mon Sep 17 00:00:00 2001
From: Martin Staab <martin.staab@aei.mpg.de>
Date: Fri, 7 Jun 2019 14:40:56 +0200
Subject: [PATCH 2/2] changed outputs of fastbinary.fourier() to long arm
outputs
---
fastsource/fastbinary/FastBinary.cpp | 16 ++++++++++++++++
fastsource/fastbinary/FastBinary.i | 17 +++++++++++------
2 files changed, 27 insertions(+), 6 deletions(-)
diff --git a/fastsource/fastbinary/FastBinary.cpp b/fastsource/fastbinary/FastBinary.cpp
index f2d1446..8f75f18 100644
--- a/fastsource/fastbinary/FastBinary.cpp
+++ b/fastsource/fastbinary/FastBinary.cpp
@@ -362,6 +362,7 @@ void FastResponse::XYZ(double f0, long q, double *XLS, double *XSL, double *YLS,
double c3 = cos(3.0*fonfs); double c2 = cos(2.0*fonfs); double c1 = cos(1.0*fonfs);
double s3 = sin(3.0*fonfs); double s2 = sin(2.0*fonfs); double s1 = sin(1.0*fonfs);
+ /*
X[2*i-1] = (c12[2*i-1] - c13[2*i-1])*c3 + (c12[2*i] - c13[2*i]) * s3 +
(c21[2*i-1] - c31[2*i-1])*c2 + (c21[2*i] - c31[2*i]) * s2 +
(c13[2*i-1] - c12[2*i-1])*c1 + (c13[2*i] - c12[2*i]) * s1 +
@@ -407,5 +408,20 @@ void FastResponse::XYZ(double f0, long q, double *XLS, double *XSL, double *YLS,
YLS[2*i] = (Y[2*i-1]*sLS + Y[2*i]*cLS);
ZLS[2*i-1] = (Z[2*i-1]*cLS - Z[2*i]*sLS);
ZLS[2*i] = (Z[2*i-1]*sLS + Z[2*i]*cLS);
+ */
+ XSL[2*i-1] = 2.0*fonfs*(c12[2*i-1]*cSL-c12[2*i]*sSL);
+ XSL[2*i] = 2.0*fonfs*(c12[2*i-1]*sSL+c12[2*i]*cSL);
+ YSL[2*i-1] = 2.0*fonfs*(c23[2*i-1]*cSL-c23[2*i]*sSL);
+ YSL[2*i] = 2.0*fonfs*(c23[2*i-1]*sSL+c23[2*i]*cSL);
+ ZSL[2*i-1] = 2.0*fonfs*(c31[2*i-1]*cSL-c31[2*i]*sSL);
+ ZSL[2*i] = 2.0*fonfs*(c31[2*i-1]*sSL+c31[2*i]*cSL);
+
+ XLS[2*i-1] = 2.0*fonfs*(c21[2*i-1]*cSL-c21[2*i]*sSL);
+ XLS[2*i] = 2.0*fonfs*(c21[2*i-1]*sSL+c21[2*i]*cSL);
+ YLS[2*i-1] = 2.0*fonfs*(c32[2*i-1]*cSL-c32[2*i]*sSL);
+ YLS[2*i] = 2.0*fonfs*(c32[2*i-1]*sSL+c32[2*i]*cSL);
+ ZLS[2*i-1] = 2.0*fonfs*(c13[2*i-1]*cSL-c13[2*i]*sSL);
+ ZLS[2*i] = 2.0*fonfs*(c13[2*i-1]*sSL+c13[2*i]*cSL);
+
}
}
diff --git a/fastsource/fastbinary/FastBinary.i b/fastsource/fastbinary/FastBinary.i
index 7b2d287..d0faf55 100644
--- a/fastsource/fastbinary/FastBinary.i
+++ b/fastsource/fastbinary/FastBinary.i
@@ -76,7 +76,8 @@ class FastGalacticBinary(lisaxml.Source):
if length == None:
length = int(T/dt)/2 + 1 # was "NFFT = int(T/dt)", and "NFFT/2+1" passed to numpy.zeros
- buf = tuple(FrequencyArray.FrequencyArray(numpy.zeros(length,dtype=numpy.complex128),kmin=kmin,df=1.0/T) for i in range(3))
+ #buf = tuple(FrequencyArray.FrequencyArray(numpy.zeros(length,dtype=numpy.complex128),kmin=kmin,df=1.0/T) for i in range(3))
+ buf = tuple(FrequencyArray.FrequencyArray(numpy.zeros(length,dtype=numpy.complex128),kmin=kmin,df=1.0/T) for i in range(6))
if status: c = countdown(len(table),10000)
for line in table:
@@ -287,12 +288,15 @@ class FastGalacticBinary(lisaxml.Source):
f0 = self.Frequency if (algorithm == 'legacy') else (self.Frequency + 0.5 * self.FrequencyDerivative * T)
if buffer == None:
- retX, retY, retZ = map(lambda a: FrequencyArray.FrequencyArray(a[::2] + 1j * a[1::2],
- dtype = numpy.complex128, # note the first frequency bin is f0 - (M/2) df,
+ #retX, retY, retZ = map(lambda a: FrequencyArray.FrequencyArray(a[::2] + 1j * a[1::2],
+ retXSL, retYSL, retZSL, retXLS, retYLS, retZLS = map(lambda a: FrequencyArray.FrequencyArray(a[::2] + 1j * a[1::2],
+ dtype = numpy.complex128, # note the first frequency bin is f0 - (M/2) df,
kmin = int(f0*T) - M/2,df = 1.0/T), # with df = 1/T
- (fastbin.XSL, fastbin.YSL, fastbin.ZSL) if (simulator == 'synthlisa') else (fastbin.XLS, fastbin.YLS, fastbin.ZLS))
+ (fastbin.XSL, fastbin.YSL, fastbin.ZSL, fastbin.XLS, fastbin.YLS, fastbin.ZLS))
+ #(fastbin.XSL, fastbin.YSL, fastbin.ZSL) if (simulator == 'synthlisa') else (fastbin.XLS, fastbin.YLS, fastbin.ZLS))
- return (retX,retY,retZ)
+ #return (retX,retY,retZ)
+ return (retXSL, retYSL, retZSL, retXLS, retYLS, retZLS)
else:
kmin, blen, alen = buffer[0].kmin, len(buffer[0]), 2*M
@@ -303,7 +307,8 @@ class FastGalacticBinary(lisaxml.Source):
# ...remember "a" is doubled up
# check: if kmin = 0, then begb = beg, endb = end, bega = 0, enda = alen
- for i,a in enumerate((fastbin.XSL, fastbin.YSL, fastbin.ZSL) if (simulator == 'synthlisa') else (fastbin.XLS, fastbin.YLS, fastbin.ZLS)):
+ #for i,a in enumerate((fastbin.XSL, fastbin.YSL, fastbin.ZSL) if (simulator == 'synthlisa') else (fastbin.XLS, fastbin.YLS, fastbin.ZLS)):
+ for i,a in enumerate((fastbin.XSL, fastbin.YSL, fastbin.ZSL, fastbin.XLS, fastbin.YLS, fastbin.ZLS)):
buffer[i][begb:endb] += a[bega:enda:2] + 1j * a[(bega+1):enda:2]
def TDI(self,T=6.2914560e7,dt=15.0,simulator='synthlisa',table=None,algorithm='mldc',oversample=1):
--
2.21.0
In order to generate time series of the galactic binary signals how the enter the long arm (science) interferometer, the lisasolve requires a small "hack" to skip the calculation of the TDI combinations X,Y,Z. Therefore, the files "FastBinary.cpp" and "FastBinary.i" in lisasolve/fastsource/fastbinary/ have to be modified! This modifiaction can be done using the `git am` command and the .patch files provided. Patch 0001 is probably not necessary but solved compatibility issues on my machine. After the modification you need to reinstall the FastBinary package.
The modification changes the definition of the `fourier()` function to return a tuple of the 6 Fourier series (H1, H2, H3, H1', H2',H3'). See example.py for more details.
import FastBinary
import numpy as np
# simulation parameters
T = 31.536e6 # one year in sec
stime = 15 # sampling time of 1 sec
mytable = np.load("/path/to/the/catalogue.npy") # for details on how to build this table see the calc_foreground directory
armlength = 2.5e9 # arm length in meters
# setup FastBinary
FastBinary.setL(armlength)
fastbin = FastBinary.FastGalacticBinary()
# run simulation and convert to time series
X = fastbin.fourier(table=mytable, T = T, dt = stime)
H1,H2,H3,h1,h2,h3 = (X[i].ifft(stime) for i in range(6))
# save output
np.save("/path/to/some/file", (H1,H2,H3,h1,h2,h3))
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment