-
Notifications
You must be signed in to change notification settings - Fork 0
/
lensing.py
executable file
·132 lines (107 loc) · 3.39 KB
/
lensing.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
import numpy as np
import time
from mpi4py import MPI
import sys
#local modules
import signalgen
import phase
import pathint
import dispath
import geopath
comm = MPI.COMM_WORLD
size = comm.Get_size()
rank = comm.Get_rank()
fmin = 311.25e6 #Hz
fban = 16e6 #Hz
nfreq = int(sys.argv[1]) # which frequency band above 311.25 to scan
freq = fmin + nfreq * fban
#fsample = 2*fban
#Period = 1.6*10**(-3) #seconds
#PulseWidth = Period/200
#R = 6.273 #c*s
fref = fmin #Hz
# conversion factors between time and points
#intime = 32/8000 # second/point
# Output prefix
FileName = "data/corrugated_sheet_"
# Generates signal
# Phase function phi(f,x)
# Path Int
# Path int scan
# Gaussian DM
# dispath
def Scan(begin, end, freq):
scan = range(begin,end)
#res = np.linspace(-1/2,1/2,3,endpoint='true')[:-1] # divides each point into subpoints of source position
res = np.array([0])
lensed = []
spec = []
for i in scan:
dp = dispath[i-3000:i+3000+1] #dispath[i-(dens*500):i+(dens*500+1)]
for j in res:
gp = geopath.generate(j)
PA = phase.PhaseArray(l,gp,dp,freq)
PI = pathint.PathInt(PA)
s1 = np.fft.irfft(sf*PI)
# save intensity as well as spectrum
lensed += [(s1**2).sum()]
spec += [sf*PI]
print(freq/10**6, i, begin, end, time.clock())
return np.array(lensed)/norm, np.array(spec)
#dispath = np.empty( 4*len(dispath.generate()) )
dispath = np.zeros(20000)
if rank == 0:
#dispath = dispath.generate()
dispath = np.load('data/corrugated_sheet.npy')
comm.Bcast(dispath, root=0)
#increasing density, do this if you're generating a new dispath
dens = 4
#temp = np.fft.rfft(dispath)
#temp = np.concatenate((temp,np.zeros( (dens-1)*2000)))
#temp = np.fft.irfft(temp)
#temp *= dens
#dispath = temp
#s = np.empty( len(Signal()) )
#if rank == 0:
# s = Signal()
#comm.Bcast(s, root=0)
#s = s[24000:27000]
# if comparing different runs, should spectrum/magnification normalize by the same signal
s = np.load('data/signal.npy')
sf = np.fft.rfft(s)
l = len(sf)
gp = geopath.generate(0)
PA = phase.PhaseArray(l,gp,gp*0,fmin)
PI = pathint.PathInt(PA)
s1 = np.fft.irfft(sf*PI)
norm = (s1**2).sum()
if rank == 0:
np.save(FileName+"Geo",gp)
np.save(FileName+"Dis",dispath)
np.save(FileName+"Signal",s)
np.save(FileName+"UnlensedSpec",sf*PI)
# number of points to scan through should be divisible by number of cores
if rank == 0:
scan = np.arange((len(gp)-1)/2, len(dispath) - (len(gp)-1)/2, (len(dispath)-(len(gp)-1)) // size )
np.save(FileName+"Scan",scan)
diff = scan[1] - scan[0]
print(scan, diff)
else:
scan = None
diff = None
scan = comm.scatter(scan, root=0)
diff = comm.bcast(diff, root=0)
scan = int(scan)
diff = int(diff)
# each processor will only process diff points for one particular frequency
mag, spec = Scan(scan,scan+diff,freq)
# comm.Barrier()
# print( "process", rank, "has", mag)
# if rank == 0 :
# magGathered = np.zeros(len(mag) * size)
# else:
# magGathered = None
#comm.Gatherv(mag, [magGathered, np.ones(size)*len(mag), np.arange(size)*len(mag), MPI.DOUBLE])
#np.save(FileName + format(freq/10**6, '.2f') + "Mag", magGathered)
np.save(FileName + format(freq/10**6, '.2f') + "Mag"+format(scan, '05') + "to" + format(scan+diff, '05'),mag)
np.save(FileName + format(freq/10**6, '.2f') + "Spec"+format(scan, '05') + "to" + format(scan+diff, '05'),spec)